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Introduction 


The  research  described  in  this  report  constitutes  a  part  of  the  long-term  scientific 
goal  of  the  Principal  Investigator,  i.e.,  to  better  understand  the  distribution  of  phyto¬ 
plankton  in  the  world’s  oceans  through  their  influence  on  the  optical  properties  of  the 
water.  Optically,  phytoplankton  reveal  their  presence  through  the  absorption  of  light  by 
their  photosynthetic  (chlorophyll  a)  and  accessory  pigments.  However,  this  absorption 
is  seldom  directly  observed  outside  a  laboratory  setting.  Instead,  it  is  usually  indirectly 
inferred  by  virtue  of  its  effect  on  the  apparent  optical  properties:  the  diffuse  reflectance 
of  the  water,  e.g.,  the  color  of  the  water,  or  the  downwelling  irradiance  attenuation  co¬ 
efficient.  Thus,  an  understanding  of  the  relationship  between  the  inherent  and  apparent 
optical  properties  is  fundamental  to  the  interpretation  of  field  observations.  Although  ul¬ 
timately  our  unders  anding  of  such  relationships  must  be  based  on  field  experiments,  the 
approach  taken  here  — The  construction  of  mathematical  models  which  simulate  as  closely 
as  possible  the  physical  processes  can  make  a  very  valuable  contribution  because  the 
various  parameters,  many  of  which  may  be  either  difficult  to  measure  or  highly  variable  in 
space  and  time,  tire  known  and  carefully  controlled  in  the  models. 

There  were  two  specific  interrelated  goals  for  the  present  research':  (1)  to  understand 
the  influence  of  the  optical  properties  of  the  ocean  as  determined  by  the  concentration  " 
of  constituents  in  Case  1  waters,  e.g.,  phytoplankton,  detrital  particles,  dissolved  organic 
material,  etc.,  on  the  transport  of  light  from  a  point  source  in  the  ocean  to  the  sea  sur¬ 
face,  for  possible  application  to  the  remote  (surface)  detection  of  bioluminescence;  and 
(2)  to  understand  the  dependence  of  the  classical  apparent  optical  properties  (the  irradi¬ 
ance  attenuation  coefficient  and  the  diffuse  reflectance)  on  the  inherent  optical  properties 
(absorption  and  scattering  coefficient  and  the  volume  scattering  function).  The  interre¬ 
lationship  between  these  goals  is  threefold:  first,  the  propagation  of  the  radiant  energy 
in  both  cases  is  governed  by  the  radiative  transfer  equation;  next,  the  same  inherent  op¬ 
tical  properties  (and  their  dependence  on  the  constituent  concentrations)  are  required  in 
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both  problems;  and  finally,  a  useful  approximate  solution  to  (1)  requires  results  from  the 
solution  to  (2). 


Summary  of  Research  Performed 

The  main  problem  of  concern  in  this  research  was  the  transport  of  radiant  energy  from 
a  point  source  in  the  ocean,  e.g.,  a  flash  of  light  emitted  by  a  bioluminescent  organism,  to 
the  surface.  The  goal  was  to  determine  a  way  of  predicting  the  horizontal  distribution  of 
irradiance  on  the  surface  from  the  spatial  distribution  and  power  of  extended  sources  in  the 
ocean  and  from  the  optical  properties  of  the  medium.  To  this  end  it  was  necessary  to  model 
the  inherent  optical  properties  of  the  water  in  terms  of  easily  measured  quantities  such 
as  the  phytoplankton  pigment  concentration,  °  as  well  as  to  solve  the  radiative  transfer 
equation  in  the  given  geometry. 

Briefly  (Gordon  1987),  a  model  of  the  optical  properties  of  the  ocean  was  developed 
providing  the  absorption  and  scattering  coefficients  of  the  medium  as  nonlinear  functions 
of  the  concentration  of  pigments  associated  with  phytoplankton  and  their  immediate  de- 
trital  material.  Monte  Carlo  computations  of  the  attenuation  coefficient  of  downwelling 
irradiance,  Kj. ,  for  an  ocean-atmosphere  system  illuminated  by  the  sun  at  zenith,  were 
found  to  agree  well  with  experimental  data,  and  thus  demonstrated  the  validity  of  the  bio- 
optical  model  for  studying  the  influence  of  phytoplankton  biomass  on  the  propagation  to 
the  surface  of  light  generated  through  binluminescence.  The  radiative  transfer  equation  for 
the  irradiance  at  the  sea  surface  resulting  .  illumination  by  a  point  source  imbedded  in 
a  homogeneous  ocean  was  then  solved  by  Monte  Carlo  techniques.  The  solution  technique 

a  By  the  term  pigment  concentration  (C)  we  mean  the  concentration  (mg/m3)  of 
chlorophyll  a  and  all  chlorophyll-like  pigments,  which  absorb  in  the  same  spectral  bands 
as  chlorophyll  o,  such  as  phaeophytin  a,  and  which  are  contained  in  phytoplankton  or  in 
their  detrital  materials.  The  sum  of  the  concentrations  of  chlorophyll  a  and  phaeophytin 
a  is  frequently  used  as  an  indicator  of  plankton  biomass. 
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was  validated  through  comparison  with  an  asymptotic  analytic  solution  for  isotropic  scat¬ 
tering.  The  computations  reveal  that  the  irradiance  distribution  just  beneath  the  surface 
as  a  function  of  R,  the  distance  measured  along  the  surface  from  a  point  vertically  above 
the  source,  is  described  by  two  regimes:  (1)  a  regime  in  which  the  irradiance  is  governed 
mostly  by  absorption  and  by  geometry,  with  scattering  playing  a  negligible  role  —  the  near 
field;  and  (2)  a  regime  in  which  the  light  field  at  the  surface  is  very  diffuse  and  the  irradi¬ 
ance  decays  approximately  exponentially  with  R  and  is  a  very  weak  function  of  the  source 
depth  —  the  diffusion  regime.  The  near  field  is  of  primary  interest  because  it  contains  most 
of  the  power  reaching  the  sea  surface.  An  analytical  model  of  the  irradiance  distribution 
just  beneath  the  surface  as  a  function  of  R,  the  source  depth,  and  the  pigment  concentra¬ 
tion  for  the  near  field  was  developed.  This  model  is  based  on  the  observation  that  at  most 
scattering  events  the  change  in  the  photon’s  direction  is  slight  and,  therefore,  scattering 
is  rather  ineffective  in  attenuating  the  irradiance.  The  analytic  solution  for  the  irradiance 
from  the  point  source  was  first  carried  out  ignoring  scattering  altogether;  however,  rec¬ 
ognizing  that  backscattering  will  attenuate  the  irradiance,  the  absorption  coefficient  was 
replaced  by  an  effective  attenuation  coefficient,  k.  This  effective  attenuation  coefficient 
was  then  determined  by  fitting  the  total  power  just  beneath  the  surface  determined  from 
the  Monte  Carlo  computations  to  the  analytical  model.  The  resulting  k  was  found  to  be 
closely  related  to  Kd ,  and  the  Monte  Carlo  irradiance  as  a  function  of  R  and  source  depth 
in  the  near  field  regime  could  be  approximated  with  high  accuracy  using  the  analytical 
model.  These  results  also  indicated  that  Kd  can  be  estimated  at  night  by  releasing  a  point 
source  in  the  water,  measuring  the  irradiance  at  the  surface  as  it  sinks,  and  fitting  the 
measurements  to  the  relationships  developed  during  the  course  of  the  research  to  deter¬ 
mine  k.  Finally,  the  analytic  model  can  also  be  inverted,  enabling  the  remote  estimation 
of  the  source  depth  and  power  from  the  irradiance  distribution  just  beneath  the  surface. 
Gordon  (1987),  provided  the  details  of  these  results,  is  reproduced  as  Appendix  1. 

The  solution  to  the  point  source  problem  above  was  incomplete  to  the  extent  that 
the  medium  was  taken  to  be  homogeneous,  i.e.,  the  optical  properties  were  independent 
of  depth.  To  remedy  this,  we  posed  the  following  question:  is  it  possible  to  estimate  the 
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point  source-generated  irradiance  in  a  stratified  ocean  given  the  vertical  distribution  of  the 
pigment  concentration  (C(z))  and  the  source  depth?  It  was  found  (Gordon  1988)  that 
the  vertical  stratification  was  unimportant  in  determining  the  horizontal  distribution  of 
irradiance  on  ihe  surface  —  the  results  of  the  homogeneous  ocean  model  could  be  extended 
to  the  stratified  ocean  if  a  constant  effective  pigment  concentration  ( C ),  defined  to  be  the 
value  of  C  for  a  homogeneous  ocean  that  would  yield  the  same  source  optical  depth  as  the 
stratified  ocean,  was  used  in  place  of  the  pigment  concentration  in  the  homogeneous  ocean 
model.  The  details  of  this  this  observation  are  presented  in  Appendix  2. 

The  point  source  problem  described  above  required  the  downwelling  irradiance  atten¬ 
uation  coefficient  just  beneath  the  surface  (Kj)  expressed  as  a  function  of  the  inherent 
optical  properties.  To  determine  this,  and  to  extend  our  understanding  of  the  relationship 
between  the  inherent  and  apparent  optical  properties  to  include  the  details  of  the  scattering 
phase  function  (volume  scattering  function  normalized  to  the  total  scattering  coefficient) 
as  well  as  environmental  factors  such  as  the  solar  zenith  angle,  the  presence  of  a  scat¬ 
tering  atmosphere,  and  sea  surface  roughness,  approximately  450  simulations  of  radiative 
transfer  in  the  ocean-atmosphere  system  were  carried  out.  Such  simulations  provide  the 
distribution  of  radiation  in  the  entire  system.  The  radiation  field  is  treated  as  experimental 
data,  albeit  data  collected  under  carefully  controlled  conditions:  a  cloud  free  sky  and  a 
homogeneous  ocean  of  precisely  known  inherent  optical  properties.  The  irradiances,  Kd , 
etc.,  can  be  derived  from  this  radiation  field  as  a  function  of  the  inherent  optical  properties 
of  the  ocean  and  ultimately  as  a  function  of  the  constituent  concentrations.  Although  the 
results  of  the  Monte  Carlo  simulations  are  still  under  analysis,  two  completed  works  have 
been  submitted  to  Limnology  and  Oceanography  for  publication. 

In  the  first  (Gordon  1989a)  it  was  shown  that  the  that  the  downwelling  irradiance 
attenuation  coefficient  just  beneath  the  surface  (AT)  and  the  mean  irradiance  attenuation 
coefficient  from  the  surface  to  the  depth  where  the  irradiance  falls  to  10%  of  its  value  at 
the  surface  (( K )  )  can  be  corrected  for  the  geometric  structure  of  the  in-water  light  field 
to  yield  quantities  that  are,  to  a  high  degree  of  accuracy,  inherent  optical  properties.  This 
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geometric  correction  is  effected  simply  by  dividing  the  experimental  irradiance  attenuation 
coefficient  by  Do,  the  downwelling  distribution  function  for  a  totally  absorbing  ocean  with 
the  same  surface  illumination,  i.e.,  Do  =  •£?od(0)/-E'tf(0))  where  i?d(0)  is  the  downwelling 
irradiance  and  £0<i(0)  is  the  downwelling  scalar  irradiance  just  beneath  the  surface.  An 
accurate  scheme  for  estimating  Do  from  simple  irradiance  and  wind  speed  measurements 
was  suggested.  The  geometry-corrected  quantities  were  then  shown  to  satisfy  the  Lambert- 
Beer  law,  i.e.,  that  Kd  can  be  linearly  resolved  into  constituent  components,  to  a  reasonable 
degree  of  accuracy,  with  the  largest  error  («  15%)  in  the  case  of  ( K )  arising  from  mixing 
nonabsorbing  particles,  e.g.,  white  sand,  with  strongly  absorbing  water  (wavelengths  > 
600  nm).  This  near-validity  of  the  Lambert-Beer  law,  when  there  are  compelling  reasons 
to  believe  that  it  should  fail,  is  shown  to  result  from  three  independent  facts:  (1)  the 
dependence  of  the  diffuse  attenuation  coefficients  on  the  geometric  structure  of  the  light 
field  can  be  removed;  (2)  pure  sea  water  is  a  much  better  absorber  than  scatterer  at 
optical  frequencies;  and  (3)  the  phase  functions  for  particles  suspended  in  the  ocean  differs 
significantly  from  that  of  pure  sea  water.  Finally,  it  is  shown  that  extrapolation  of  the 
corrected  diffuse  attenuation  coefficients  to  the  limit  c  — +  cw,  where  c  is  the  total  beam 
attenuation  coefficient  and  cw  is  the  beam  attenuation  coefficient  of  pure  sea  water,  yields 
quantities  that  are  within  2%  of  the  corresponding  quantities  that  would  be  measured 
for  an  ocean  consisting  of  pure  sea  water  with  the  sun  at  the  zenith  and  the  atmosphere 
removed.  The  details  of  the  analysis  of  Kd  are  presented  in  Appendix  3. 

In  the  second  (Gordon  1989b,  and  Appendix  4),  the  variation  of  the  diffuse  reflectance 
of  natural  waters  with  sun  angle  was  studied  is  found  to  be  dependent  on  the  shape 
of  the  volume  scattering  function  (VSF)  of  the  medium.  It  was  also  found  that  single 
scattering  theory  can  be  used  to  estimate  the  reflectance  -  sun  angle  variation  given  the 
VSF,  and  conversely,  the  VSF  can  be  retrieved  from  measurements  of  the  variation  of 
the  reflectance  with  sun  angle.  The  complex  variation  of  reflectance  with  the  incident 
illumination  and  surface  roughness  was  found  to  be  reducible  to  the  variation  of  a  single 
parameter:  R(D0)  =  kD0R(l),  where  A:  is  a  constant,  and  R(  1)  is  the  reflectance  the  ocean 
would  have  with  the  atmosphere  removed  and  the  sun  at  the  zenith,  i.e.,  the  variation  of 
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R  with  the  incident  illumination  and  the  surface  roughness  can  be  completely  explained 
through  their  effect  on  D0,  which  was  defined  in  the  preceeding  paragraph.  The  parameter 
k  depends  mostly  on  the  scattering  phase  function,  and  for  Wo  <  0.9,  where  u>o  is  the  ratio 
of  the  scattering  coefficient  to  the  beam  attenuation  coefficient  of  the  medium,  it  can  be 
computed  using  the  single  scattering  approximation.  These  observations  are  applicable  to 
all  but  the  most  reflective  of  natural  waters. 

In  addition,  two  other  papers  were  published  (Gordon  and  Castaiio  1988,  Gordon 
and  Castano  1989)  which,  although  they  had  no  bearing  on  the  specific  problems  under 
consideration  here,  were  completed  while  Diego  Castano  was  supported  as  a  Graduate 
Assistant  on  the  Contract  and  thus  acknowledgement  of  ONR  support  for  these  papers 
is  appropriate.  The  first  paper  shows  that  the  atmospheric  correction  algorithm  for  the 
Coastal  Zone  Color  Scanner  (CZCS)  should  be  largely  unaffected  by  the  eruption  of  the 
volcano  El  Chichon  in  late  March  -  early  April  1982,  while  the  second  provides  a  novel 
method  for  examining  oceanic  aerosols  using  CZCS  imagery.  For  completeness,  these  are 
reproduced  in  Appendices  5  and  6  respectively. 
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the  surface.  The  goal  was  to  determine  a  way  of  predicting  the  horizontal  distribution  of 
irradiance  on  the  surface  from  the  spatial  distribution  and  power  of  extended  sources  in  the 
ocean  and  from  the  optical  properties  of  the  medium.  To  this  end  it  was  necessary  to  model 
the  inherent  optical  properties  of  the  water  in  terms  of  easily  measured  quantities  such 
as  the  phytoplankton  pigment  concentration,  0  as  well  as  to  solve  the  radiative  transfer 
equation  in  the  given  geometry. 

Briefly  (Gordon  1987),  a  model  of  the  optical  properties  of  the  ocean  was  developed 
providing  the  absorption  and  scattering  coefficients  of  the  medium  as  nonlinear  functions 
of  the  concentration  of  pigments  associated  with  phytoplankton  and  their  immediate  de- 
trital  material.  Monte  Carlo  computations  of  the  attenuation  coefficient  of  downwelling 
irradiance,  Kj,  for  an  ocean-atmosphere  system  illuminated  by  the  sun  at  zenith,  were 
found  to  agree  well  with  experimental  data,  and  thus  demonstrated  the  validity  of  the  bio- 
optical  model  for  studying  the  influence  of  phytoplankton  biomass  on  the  propagation  to 
the  surface  of  light  generated  through  bioluminescence.  The  radiative  transfer  equation  for 
the  irradiance  at  the  sea  surface  resulting  from  illumination  by  a  point  source  imbedded  in 
a  homogeneous  ocean  was  then  solved  by  Monte  Carlo  techniques.  The  solution  technique 

a  By  the  term  pigment  concentration  (C)  we  mean  the  concentration  (mg/m3)  of 
chlorophyll  a  and  all  chlorophyll-like  pigments,  which  absorb  in  the  same  spectral  bands 
as  chlorophyll  a,  such  as  phaeophytin  a,  and  which  are  contained  in  phytoplankton  or  in 
their  detrital  materials.  The  sum  of  the  concentrations  of  chlorophyll  a  and  phaeophytin 
a  is  frequently  used  as  an  indicator  of  plankton  biomass. 
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Bio-optical  model  describing  the  distribution  of  irradiance 
at  the  sea  surface  resulting  from  a  point  source 
embedded  in  the  ocean 


Howard  R.  Gordon 


A  model  of  the  optical  properties  of  the  ocean,  providing  the  absorption  and  scatterh-g  Loefrcienta  of  the 
medium  as  nonlinear  functions  of  the  concentration  of  pigments  associated  with  phytoplankton  and  their 
immediate  detrital  material,  is  presented.  Monte  Carlo  computations  of  the  attenuation  coefficient  of 
downwelling  irradiance  Ka  for  an  ocean-atmosphere  system  illuminated  by  the  sun  at  zenith,  agree  well  with 
experimental  data  and  demonstrate  the  validity  of  such  a  model  for  studying  the  influence  of  phytoplankton 
biomass  on  the  propagation  to  the  surface  of  light  generated  through  bioluminescence.  The  radiative 
transfer  equation  for  the  irradiance  at  the  sea  surface  resulting  from  illumination  by  a  point  source  embedded 
in  the  water  is  solved  by  Monte  Carlo  techniques.  The  solution  technique  is  validated  through  comparison 
with  an  asymptotic  analytic  solution  for  isotropic  scattering.  The  computations  show  that  the  irradiance 
distribution  just  beneath  the  surface  as  a  function  of  R,  the  distance  measured  along  the  surface  from  a  point 
vertically  above  the  source,  is  described  by  two  regimes:  (1)  a  regime  in  which  the  irradiance  is  governed 
mostly  by  absorption  and  geometry  with  scattering  playing  a  negligible  role — the  near  field;  (2)  a  regime  in 
which  the  light  field  at  the  surface  is  very  diffuse  and  the  irradiance  decays  approximately  exponentially  in  R 
and  is  a  very  weak  function  of  the  source  depth — the  diffusion  regime.  The  near  field  is  of  primary  interest 
because  it  contains  most  of  the  power  reaching  the  sea  surface.  An  analytical  model  of  the  irradiance 
distribution  just  beneath  the  surface  as  a  function  of  R,  the  source  depth,  and  the  pigment  concentration  for 
the  near  field  is  presented.  This  model  to  based  on  the  observation  that  at  most  scattering  events  the  change 
in  the  photon’s  direction  is  slight,  and  therefore,  scattering  to  rather  ineffective  in  attenuating  the  irradiance. 
An  analytic  solution  for  the  irradiance  from  the  point  source,  then,  to  first  carried  out  ignoring  scattering 
altogether;  however,  recognizing  that  backscattering  will  attenuate  the  irradiance,  the  absorption  coefficient 
is  replaced  by  an  effective  attenuation  coefficient  k.  This  effective  attenuation  coefficient  to  determined  by 
fitting  the  total  power  just  beneath  the  surface  determined  from  the  Monte  Carlo  computations  to  the 
analytical  model.  The  resulting  k  to  closely  related  to  Kj,  and  the  Monte  Carlo  irradiance  as  a  function  of  R 
and  source  depth  in  the  near-field  regime  can  be  approximated  with  high  accuracy  using  the  model.  These 
results  indicate  Kj  can  be  estimated  at  night  by  releasing  a  point  source  in  the  water,  measuring  the  irradiance 
at  the  surface  as  it  sinks,  and  fitting  the  measurements  to  the  relationships  developed  here  to  determine  k. 
The  analytic  model  also  enables  estimation  of  the  source  depth  and  power  from  the  irradiance  distribution 
just  beneath  the  surface. 


I.  Introduction 

In  an  effort  to  examine  the  extent  to  which  marine 
bioluminescent  emissions  can  be  studied  from  ships 
and/or  aircraft  it  is  necessary  to  understand  the  sur¬ 
face  manifestation  of  the  bioluminescent  signal  in 
terms  of  the  optical  properties  of  the  water  and  the 
strength  of  the  emission.  Thus  the  problem  examined 
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here  consists  of  determining  the  distribution  of  irradi¬ 
ance  at  the  sea  surface  due  to  an  isotropically  emitting 
point  source  embedded  in  the  ocean  and  the  depen¬ 
dence  of  this  irradiance  on  the  optical  properties  of  the 
water,  i.e.,  given  a  distribution  of  isotropically  emitting 
point  sources  embedded  in  the  ocean,  the  variation  of 
the  associated  light  field  at  the  sea  surface  with  the 
optical  properties  of  the  water  is  determined.  Since 
the  duration  of  bioluminescent  flashes  tends  to  be  long 
compared  to  the  mean  lifetime  of  a  photon  in  the  water 
(absorption  mean  free  path/speed  of  light)  of  at  most 
~1  na,  the  light  field  can  be  described  by  steady-state 
radiative  transfer  theory.  The  steady-state  radiative 
transfer  equation  (RTE),  which  describes  the  trans¬ 
port  of  the  radiance  [L(r,{ )]  at  r  in  a  direction  specified 
by  the  unit  vector  (,  is  given  by 
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({■  V)/.(ri)  +  dr)Ur,{) 

=  [  /»b\f'  -  ()Ur,(')dU(?)  +  Q(r,{),  (X) 
hr 

where  /9(r,£'  -*  £)  is  the  volume  scattering  function  at  r 
for  scattering  from  the  direction  £'  into  the  direction  |, 
c(r)  is  the  beam  attenuation  coefficient  at  r,  and  Q(r,£) 
j  is  the  intensity  (power  per  unit  solid  angle)  per  unit 
volume  of  sources  at  r  in  the  direction  |.  c(r)  can  be 
related  to  the  absorption  coefficient  a( r)  and  the  scat¬ 
tering  coefficient  b(r)  by 

c(r)  =  a(r)  +  6(r),  (2) 

where 

6(r)  «=  [  0(r,('  -  Odntf’h  (3) 

Jtr 

It  can  be  shown1  that  for  a  closed  volume  V  with  «o(r) a 
:  b(r)/c(r)  <  1,  surrounded  by  a  surface  S,  the  RTE 

possesses  unique  solutions  given  the  sources  Q(r,() 
within  V and  the  radiance  into  V  from  the  outside,  i.e., 
[L(rs,|)J  for  £  •  h  <  1  (defined  to  be  (L(inc,(rs>|)j),  where 
r,  is  on  the  bounding  surface  S,  and  h  is  the  outward 
unit  normal  to  S.  In  the  problem  of  interest  here,  the 
energy  source  is  isotropically  emitting  radiation  of  unit 
intensity  at  a  point  r0  in  the  water,  and  there  is  no 
energy  incident  from  the  sea  surface,  i.e., 

Q(r,|)  =  «(r  -  r0), 

L,inet(rJ)  =  0. 


Thus  the  quantities  a(r)  and  /9(r,£'  —  |)  are  all  that  are 
required  to  predict  the  transport  of  the  radiant  energy 
to  the  sea  surface.^  Traditionally,  the  volume  scatter¬ 
ing  function  -*  |)  is  normalized  to  the  local 

scattering  coefficient,  replacing  it  by  the  scattering 
phase  function  P(r,£'  -»  £),  defined  according  to 

P( r.t  -  I)  =  0(r,f'  -  i)tb(T). 

Assuming  that  the  optical  properties  of  the  medium 
are  independent  of  position  and  introducing  these  def¬ 
initions  into  the  RTE,  we  have 

({ •  v)Hr£)  +  cL(tJ)  ; 

=  b  f  P(i’  -  j)L(r.l')dnd')  +  a(r  -  r0).  (4) 

» 

Dividing  by  c  and  recalling  that  5(ax)  =  a-1$(x)  yields 

U •  V')L(cr,|)  +  L(cr,{)  j 

”  “o  f  f)i(cr ,(')dn((')  +  S(cr  -  cr0),  (5) 

h 

where  the  prime  on  the  gradient  operator  indicates 
that  the  derivatives  ape  now  with  respect  to  the  optical 
variables  cx,  cy,  and  at.  Using  the  optical  variables  we 
see  that  the  radiance  in  the  medium  is  governed  com¬ 
pletely  by  «0  and  P( £'  -*■  £). 

It  is  important  to  note  that  the  scattering  and  ab¬ 
sorption  properties  of  the  medium  can  be  conveniently 
separated  into  those  due  to  the  water  itself  and  those 
due  to  materials  dissolved  or  suspended  in  water,  i.e.. 


0  “  a»  +  X  a“ 

I 

4 

and  so 


b  =  bu  + 


bP  -bPw  +  Y  bP„ 

I 

where  the  subscript  w  refers  to  the  water  itself,  and  the 
subscript  i  refers  to  the  ith  constituent  of  the  medium. 
It  is  necessary  then  to  determine  t»o  and  P,  given  aw,  o„ 
du,  and  di*. 

We  will  consider  only  a  two-component  system  con¬ 
sisting  of  water  and  phytoplankton  (and  their  immedi¬ 
ate  derivatives),  i.e., 


o  =  ow  +  ap 

+  or  bP  -  bwPw  + 


where  the  subscript  p  refers  to  phytoplankton.  Given 
the  above  quantities,  w0  and  P  can  be  determined  from 


ujcjcj  + 

w°m  -  c/tr+T  • 

u pPpkp/cJ  +  U>UPW 
WaP — — 


The  choice  of  these  parameters  and  the  particle  scat¬ 
tering  phase  function  is  discussed  in  the  next  section. 


H.  Optical  Model  of  Water  and  Its  Constituents 

To  provide  a  realistic  computation  of  the  transport 
of  radiant  energy  from  a  point  source  in  the  ocean,  a 
model  for  the  optical  properties  of  the  medium  is  need¬ 
ed.  Such  a  model  must  provide  the  scattering  phase 
and  absorption  coefficients  of  the  water  and  its  con¬ 
stituents.  In  waters  labeled  Case  1  by  Morel  and 
Prieur2  (and  comprising  much  of  the  open  ocean)  the 
optical  properties  are  controlled  by  phytoplankton 
and  their  associated  detrital  material.  Thus  the  mod¬ 
el  should  provide  the  optical  properties  of  the  phyto¬ 
plankton  and  their  associated  detrital  material  as  a 
function  of  some  convenient  measure  of  the  biological 
activity  in  the  water.  The  usual  indicator  of  the  bio¬ 
logical  activity  is  the  concentration  of  chlorophyll  a.3 
This  is  taken  to  be  a  measure  (although  imprecise)  of 
the  phytoplankton  biomass.  In  Case  1  waters  it  is  to 
be  expected,  that  to  the  extent  the  phytoplankton 
concentration  can  be  parametrized  by  the  chlorophyll 
a  concentration,  the  optical  properties  of  the  biological 
material  can  also  be  specified  by  the  chlorophyll  a 
concentration.  This  has  been  found  to  be  the  case, 
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although  this  specification  can  only  be  made  in  a  sta¬ 
tistical  sense  due  to  the  natural  variability  of  the  rela¬ 
tive  concentrations  of  the  absorbing  pigments  associ¬ 
ated  with  phytoplankton. 


A.  Optical  Properties  of  Pure  Seawater 

The  absorption  coefficient  of  pure  seawater  has  been 
inferred  from  measurements  of  downwelling  and  up- 
welling  irradiance  in  oligotrophic  waters  such  as  the 
Sargasso  Sea.2  4  The  scattering  coefficient  has  been 
measured  directly  for  pure  water  and  for  saline  solu¬ 
tions  of  pure  water  corresponding  to  salinities  between 
35  and  39  ppt  by  Morel.5  The  results  of  these  mea¬ 
surements  are  presented  in  Table  I.  The  scattering 
phase  function  for  pure  water  is  given  by  the  familiar 
Rayleigh  scattering  formula: 


Morel’s  measurements  suggest  a  depolarization  factor 
5  of  ~0.09;  thus 

Pw(e)  =  Pu(90°)(l  +  0.835  COST0).  (9) 


B.  Optical  Properties  of  the  Particles 

The  scattering  coefficient  of  particles  in  the  ocean 
has  been  studied  as  a  function  of  the  pigment  concen¬ 
tration  by  Morel6  (also  see  Ref.  7.)  The  result  of 
measurements  in  Case  1  waters  from  several  locations 
indicates  that  the  scattering  coefficient  at  550  nm, 
M550),  is  nonlinearly  related  to  the  pigment  concen¬ 
tration  C  through 

bc  =  BcC°62,  (10) 

where  6c(550)  is  in  m-1  and  C  is  in  mg/m3.  The 
constant  Be,  the  scattering  coefficient  at  a  pigment 
concentration  of  1  mg/m3,  varies  from  0.12  to  0.45  with 
an  average  value  of  0.30.  This  variation  in  Be  is  due  to 
the  natural  variability  of  scattering  over  the  various 
species  of  phytoplankton  as  well  as  a  variability  in 
scattering  by  the  detrital  particles  associated  with  the 
phytoplankton. 

The  absorption  of  phytoplankton  and  their  associat¬ 
ed  detrital  material  (excluding  yellow  substances)  has 
been  deduced  from  irradiance  measurements  by 
Prieur  and  Sathyendranath.®  Their  study  indicated  a 
nonlinear  relationship  between  absorption  and  pig¬ 
ment  concentration  similar  to  that  described  by  Smith 
and  Baker.910  The  source  of  the  nonlinearity  in  the 
absorption-pigment  concentration  relationship  (and 
also  the  scattering-pigment  concentration  relation¬ 
ship)  is  believed  to  be  a  systematic  variation  in  the 
ratio  of  the  concentration  of  phytoplankton  to  that  of 
detrital  material  as  a  function  of  tbe  concentration  of 
phytoplankton.  Hobson  et  al.u  observed  that  as  the 
concentration  ,of  phytoplankton  increases  (and  thus 
the  pigment  concentration  also)  the  ratio  of  phyto¬ 
plankton  carbon  to  detrital  carbon  also  increases. 
Thus,  in  the  two-component  absorption  system  of  phy¬ 
toplankton  and  their  detrital  material,  the  relative 
amounts  of  the  components  vary  with  the  pigment 
concentration  forcing  the  total  absorption  to  be  a  non- 


TaU*  I.  Abaorption  and  Scattering  Coaftlclanta  of  Puia  Saawatar 


X 

(nm) 

K 

<m-') 

o„. 

(m-1) 

U)u, 

420 

0.0061 

0.0153 

0.285 

440 

0.0049 

0.0145 

0.253 

460 

0.0041 

0.0156 

0.208 

480 

0.0034 

00176 

0.162 

500 

0.0029 

0.0257 

0.101 

520 

0.0024 

0.0477 

0.048 

540 

0.0021 

0.0558 

0.036 

550 

0.0019 

0.0638 

0.029 

560 

0.0018 

0.0708 

0.025 

580 

0.0016 

0.108 

0.015 

600 

0.0014 

0.244 

0.006 

620 

0.0012 

0.309 

0.004 

linear  function  of  the  pigment  concentration.  This 
nonlinearity  can  be  molded  in  many  ways.  Smith  and 
Baker9-10  used  a  linear  relationship  between  absorp¬ 
tion  (diffuse  attentuation  coefficient  Kd)  and  pigment 
concentration  but  with  different  slopes  ( dKJdC )  above 
and  below  1  mg/m3.  Prieur  and  Sathyendranath® 
found  that  a  similar  relationship  could  be  used  to 
explain  their  absorption  data  but  also  discovered  that 
for  C  <  10  mg/m3  a  single  equation 

ac(X)  =0.06^c(X)C0602,  (11) 


where  ac(X)  is  in  m_1  and  C  is  in  mg/m3,  fit  the  experi¬ 
mental  data  as  well  as  the  segmented  linear  relation¬ 
ship.  In  this  equation  Ac(X)  is  the  absorption  coeffi¬ 
cient  of  phytoplankton  normalized  to  440  nm,  i.e., 


AC(X)  = 


ac(M 
M440)  ’ 


The  relative  absorption  of  phytoplankton  Ac(X)  de¬ 
duced  by  Prieur  and  Sathyendranath  agrees  well  with 
absorption  measurements  made  on  phytoplankton 
cultures  by  Sathyendranath.12  Note  that  ac(X)  is  the 
absorption  coefficient  of  “phytoplankton  and  their  im¬ 
mediate  derivatives  or  by-products  having  similar  op¬ 
tical  effects.”®  We  assume  here  that  there  are  no  other 
sources  of  particles  in  the  water,  such  as  resuspended 
particles  from  the  bottom  in  coastal  areas;  i.e.,  we  limit 
the  discussion  to  Case  1  waters,  and  as  such  ac(X) 
represents  the  absorption  coefficient  of  all  the  parti¬ 
cles  in  the  water.  A  similar  meaning  is  also  attached  to 
MX). 

It  is  interesting  to  note  that  MX)  and  ac(X)  vary 
with  pigment  concentration  in  nearly  the  same  man¬ 
ner,  i.e.,  approximately  as  C06.  This  suggests  that 
MX)/ac(X)  is  nearly  independent  of  the  pigment  con¬ 
centration,  and  in  fact 


^c(M  SfiM 

- *  1 6.6 - • 

oc(X)  AC(X) 


(12) 


This  provides  an  estimation  of  wp(X)  and  shows  that 
this  quantity  is  in  the  first  approximation  independent 
of  the  pigment  concentration.  Table  II  gives  o>p(X)  at 
440  and  550  nm  for  the  three  different  values  of 
Bc(550):  the  minimum,  mean,  and  maximum  ob¬ 
served  in  Case  1  waters.  To  estimate  «p(440)  it  is 
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T»W*  II.  «p(\)  (or  Seawater  P article*  Auumlng  Bc(A)  =  Bc<5$0)  X 
(550/A)" 


Bt  (550) 

wp(440) 

wp(550) 

n  =  -1 

n  -  0 

n  =  1 

0.12 

0.606 

0.660 

0.713 

0.848 

0.30 

0.800 

0.833 

0.862 

0.933 

0.45 

0.903 

0.954 

necessary  to  assume  a  wavelength  dependence  for 
Bc(X).  The  computations  presented  in  Table  II  of 
oip(440)  have  been  made  assuming  Bc(\)  ~  X_n  with  n 
=  -1,  0,  and  1.  Since  the  scattering  by  absorbing 
particles  tends  to  be  depressed  in  the  absorption 
bands,13  values  of  n  =  —1  to  0  will  be  favored  for 
phytoplankton  in  the  440-550-nm  spectral  region.14 
Table  III  provides  values  of  wp(X)  for  four  species  of 
cultured  phytoplankton  studied  by  Bricaud  et  al.13 
Clearly,  values  of  Bc(550)  can  be  found  which  will 
bring  o>p(X)  into  conformity  with  the  culture  measure¬ 
ments.  One  must  remember,  however,  that  i>c(X)  and 
qc(X)  include  the  contribution  of  any  detritus  which 
covaries  with  C.  The  apparent  lack  of  any  dependence 
of  c op  on  C  suggests  that  the  plankton  and  their  imme¬ 
diate  detritus  must  have  very  similar  absorption  and 
scattering  properties  at  least  in  a  statistical  sense. 

To  model  the  transport  of  bioluminescent  emission 
in  the  water,  the  optical  properties  of  the  medium  near 
the  wavelengtn  of  maximum  emission  are  required. 
This  is  centered  near  480  nm.  The  data  in  Table  II  for 
the  mean  value  of  J3C(550)  indicate  that  up  is  likely  to 
fall  in  the  0.80-0.94  range.  The  actual  choice  «p  re¬ 
quires  a  choice  for  n;  however,  the  value  of  «p  is  very 
insensitive  to  n,  ranging  from  0.845  to  0.891  as  n  varies 
from  - 1  to  2.  Since,  for  the  most  part,  we  are  interest¬ 
ed  in  pigment  concentrations  smaller  than  1  mg/m3, 
for  which  a  significant  portion  of  the  scattering  will  be 
due  to  the  detrital  material,  we  choose  n  =  1,  yielding  a 
nominal  value  of  0.88  for  «p(480).  Also,  we  have  cho¬ 
sen  0.16,  the  value  at  480  nm,  for  ww. 

The  remaining  parameters  needed  to  specify  the 
optical  properties  of  the  medium  are  cp/cw  and  the 
particle  phase  function.  From  Table  I  it  is  seen  that  cw 
-  0.021  m-1  at  480  nm.  Noting  that  Ac(480)  =  0.798 
and  using  the  mean  value  of  Bc(5 50)  with  a  X-1  spectral 
variation,  we  have 

c„  «  (0.34  +  0.048)C°e, 


*  18.5C06,  (13) 

Cu> 

at  480  nm.  This  provides  the  connection  between  cp/ 
cw  and  the  pigment  concentration. 

The  particle  phase  function  is  the  most  difficult 
quantity  to  parametrize,  because  it  requires  the  indi¬ 
vidual  phase  functions  of  the  plankton  and  the  detrital 
material  and  their  respective  scattering  coefficients. 
This  presents  a  serious  problem,  since  neither  phase 
function  has  ever  been  measured.  Thus  we  are  forced 
to  rely  on  measurements  of  the  total  particle  function 


Table  III.  MjX)  lot  Four  Sp*dw  ot  Phytoplankton  (altar  Bricaud  #1  at’*) 


Species 

taip(440) 

o»p(550) 

Hymenomonas  elongate 

0.65 

0.89 

Platymonas  Sp 

0.76 

0.90 

Tetraselmis  maculate 

0.76 

0.91 

Coccolithua  huxleyi 

0.88 

0.97 

SCATTERING  ANGLE  6  (Deg.) 

Fig.  1.  Particle  scattering  phase  function  used  in  this  study.  The 
water  phase  function  (lower  curve  at  small  scattering  angles)  is  also 
shown  for  comparison. 


(plankton  plus  detrital  material).  Petzold15  has  mea¬ 
sured  volume  scattering  functions  at  530  nm  for  waters 
in  several  locations  with  very  different  turbidities  (to¬ 
tal  scattering  coefficients).  After  the  scattering  by 
pure  water  is  subtracted,  the  resulting  particle  phase 
functions  have  a  standard  deviation  which  is  within 
~30%  of  the  mean.  This  is  remarkable  considering 
that  the  particle  scattering  coefficient  varied  over  a 
factor  of  ~50.  The  mean  particle  phase  function  and 
its  standard  deviation  are  shown  in  Fig.  1  along  with 
the  phase  function  for  scattering  by  the  water  itself 
lPu,(6)].  This  mean  particle  phase  function  derived 
from  Petzold’s  measurements  is  adopted  for  this 
study.  It  should  be  pointed  out,  however,  that  this 
phase  function  will  be  adequate  only  for  simulations 
such  as  those  under  investigation  in  this  study,  which 
are  expected  to  be  insensitive  to  the  particle  backscat- 
tering  probability  (6b)p  defined  according  to 

(6»)p  E  (bh)p/bp,  (14a) 

where 

bb  =  2r  (  0(0)  sinOdQ.  (14b) 

This  restriction  arises  because  (6b)p  measured  for  oli- 
gotrophic  waters,  where  the  concentration  of  detrital 
material  is  relatively  high  compared  with  phytoplank¬ 
ton,  is  significantly  higher  than  (Bb)p  measured  in 
euthrophic  waters,  where  the  detrital  concentration  is 
low,  and  in  plankton  cultures  where  the  detritus  is 
essentially  absent.  Thus  a  single  phase  function,  inde¬ 
pendent  of  the  viable-phytrtpla-'kton-tc-dctrital-par- 
ticle  ratio,  cannot  be  expected  to  be  useful  for  compu- 
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tations  of  the  diffuse  reflectance  of  the  ocean,  which  is 
proportional  to  bb/a.1 

C.  Range  of  Validity  of  the  Optical  Model 

To  provide  a  measure  of  the  extent  to  which  this 
optical  model,  relating  the  inherent  optical  properties 
of  the  ocean  to  the  pigment  concentration,  can  repro¬ 
duce  the  apparent  optical  properties  of  the  ocean,  the 
RTE  has  been  solved  for  sunlight  at  480  nm  incident  on 
the  top  of  the  atmosphere  from  the  zenith  to  obtain  the 
downwelling  irradiance  attenuation  coefficient  K<j(z) 
defined  by 


where  Ed(z)  is  the  downwelling  irradiance  at  a  depth  z. 
Rayleigh  and  aerosol  scattering  are  both  accounted  for 
in  the  simulation;  however,  the  sky  is  assumed  to  be 
cloud  free.  The  cotaputations  were  carried  out  using 
Monte  Carlo  techniques  for  pigment  concentrations 
varying  from  0  to  about  4.5  mg/m3.  The  computations 
at  C  =  0  yield  the  diffuse  attenuation  coefficient  for 
pure  seawater  Kw  at  480  nm,  which  was  determined  to 
be  0.0206  m“l,  in  reasonably  good  agreement  with  the 
0.0194  m-1  estimated  by  Baker  and  Smith.4  After  the 
coefficient  for  pure  Water  was  subtracted,  an  excellent 
fit  to  the  data  for  the  mean  Kd  —  Kw  over  a  depth 
interval  of  about  half  of  the  euphotic  depth  (z  =  0  to  z  *• 
2.3 lKd)  to 

Kd  —  Kw  =  0.070C0615  (16a) 

was  obtained.  This  should  be  compared  with  Morel’s 
(see  Ref.  7)  fit  of  the  same  quantity  to  actual  experi¬ 
mental  data,  which  yielded 

Kd-Kw  =  0.074 C0703  (16b) 

For  the  range  in  pigment  concentration  between  0.05 
and  2  mg/m3  these  two  results  agree  within  20%,  allow¬ 
ing  the  determination  of  Kd  to  better  than  5%  for  0  ^  C 
<  1  mg/m3  and  better  than  10%  for  0  <  C  <  2  mg/m3. 
The  comparison  between  the  measurements  and  the 
model  is  much  poorer  for  larger  C,  reaching  ~30%  for  C 
near  10  mg/m3.  This  agreement  between  the  present 
model  and  the  Morel  measurements  should  not  be 
surprising,  considering  that  part  of  the  model  itself  is 
based  on  an  empirical  relationship  between  a  and  C 
which  was  derived  .from  the  Morel  data  set.8  The 
agreement  however,  does  attest  to  the  correctness  of 
the  empirical  6  —  C  relationship  and  the  reasonable¬ 
ness  of  the  scattering  phase  function  used  in  the  com¬ 
putations. 

The  standard  approximation  Kd  =*  a  +  ftj,  yielded 

Kd-Kw^  0.052C06,  (16c) 

while  the  Monte  Carlo  computations  yielded 

Kd-Kw- 0.057 C06'7,  (16d) 

when  the  irradiance  attenuation  coefficients  were  eval¬ 
uated  just  beneath  the  surface.  The  agreement  be¬ 
tween  the  approximation  and  the  exact  computations 
near  the  surface  attests  to  the  accuracy  of  the  approxi¬ 
mation.  It  must  be  noted,  however,  that  Kd  is  a  func- 


Fig.  2.  Physical  setting  for  defining/?,  z0, 9,  and£.  The  sea  surface 
is  in  the  x-y  plane.  Note  that  the  light  field  in  the  medium  is 
independent  of  4>  for  an  isotropic  point  source  at  zq. 

tion  of  depth,  while  the  standard  approximation  re¬ 
sults  in  a  constant  Kd.  The  basic  difference  between 
the  surface  value  of  Kd  and  that  measured  over  a  depth 
range  z  =  0  to  z  »  2.3 /Kd  is  an  increase  of  the  value  of 
Kd  —  Kw  at  C  =  1  mg/m3;  the  functional  dependence  on 
C  is  essentially  unchanged. 

The  ability  of  the  model  to  reproduce  the  Kd  —  C 
relationship  derived  from  Morel's  data7  suggests  that 
the  model  of  the  inherent  optical  properties  as  a  func¬ 
tion  of  the  pigment  concentration  will  be  useful  in 
studying  the  propagation  of  light  emitted  from  a  point 
source  in  the  ocean. 

HI.  Solution  of  the  RTE  for  a  Point  Source  In  the  Ocean 

Most  standard  techniques  for  solving  the  RTE  re¬ 
quire  that  the  spatial  dependence  of  the  light  field  be 
described  by  a  single  coordinate,  i.e.,  in  the  case  of  the 
ocean  illuminated  by  solar  radiation,  the  radiance  is 
assumed  to  be  a  function  of  depth  only,  not  the  hori¬ 
zontal  coordinates.  Clearly,  in  the  case  of  point- 
source  illumination,  the  radiance  depends  on  both 
depth  and  horizontal  position  [Eq.  (4)j,  and  the  stan¬ 
dard  solution  techniques  will  no  longer  be  applicable. 
Thus  in  this  study  the  RTE  is  solved  using  Monte 
Carlo  techniques,  which  are  applicable  to  all  geome¬ 
tries  and  which  can  easily  be  modified  to  include  both 
vertical  and  horizontal  variations  of  the  optical  prop¬ 
erties  of  the  water  if  desired. 

A.  Monte  Carlo  Technique 

Given  an  isotropically  emitting  point  source  in  the 
ocean,  we  wish  to  calculate  the  spatial  and  angular 
distribution  of  the  irradiance  transmitted  to  the  sea 
surface  as  a  function  of  the  pigment  concentration  in 
the  water.  Since  all  the  inherent  optical  properties  of 
the  medium  have  been  given  above  as  a  function  of  C,  it 
is  only  necessary  to  solve  the  RTE  to  determine  the 
desired  distributions.  The  solution  is  effected  by 
Monte  Carlo  techniques.  Photons  are  emitted  isotro¬ 
pically  from  a  point  at  depth  z0  beneath  the  surface 
(see  Fig.  2).  This  emission  can  be  accomplished  by 
choosing  the  polar  angle  0  (the  angle  between  the  pho¬ 
ton’s  direction  and  the  nadir)  from  cos0  =  1  —  pu  where 
Pi  is  a  random  number  distributed  uniformly  on  [0,1], 
and  choosing  the  azimuth  angle  4>  according  to  4>  *= 
2*-p2,  where  P2  is  a  second  random  number  with  the 
same  distribution  as  pi.  However,  in  the  ocean,  scat¬ 
tering  is  mostly  in  the  near-forward  direction  so  pho- 
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tons  emitted  downward  from  the  source,  i.e.,  away 
from  the  interface,  have  little  chance  of  reaching  the 
surface  to  contribute  to  the  irradiance.  Following 
these  photons  for  all  practical  purposes,  a  waste  of 
computational  time.  Therefore,  in  the  present  code 
we  have  chosen  to  sample  6  from  a  distribution  which 
yields  an  increased  number  of  photons  starting  toward 
the  surface.  The  normalized  probability  density  for 
this  distribution  is 


pf(o )  = 


Ji~-  <2'  , 

jt(  1  +  f  cosfl) 


(17) 


where  0  <  <  <  1,  and  the  subscript  /  indicates  false  in 
that  the  true  probability  density  from  which  8  should 
be  chosen  is  pt(6 )  =  */>  sintf.  Sampling  8  from  the 
incorrect  probability  distribution  introduces  a  bias 
which  can  be  removed  by  providing  each  photon  with 
an  initial  statistical  weight  W  given  by 


W 


p,{H) 

P,cn 


Using  p/(8)  above,  given  a  random  number  pu  the  angle 
9  is  chosen  from 

l.  V  lot.  -  — “  •  (18) 

1  +  V- 


where 


>•  =* 


(19) 


and  the  weight  of  the  associated  photon  is  initialized  to 


1  +  ( cusu  .  _ 

IV  =  x  -  -  —  — -  sinS  •  (20) 

2\l  -  t2 

The  photon’s  path  is  then  followed  using  standard 
Monte  Carlo  techniques.  At  each  collision  its  weight 
is  multiplied  wo-  On  encountering  the  surface  its  hori¬ 
zontal  position  R  is  determined,  and  the  associated 
upwelling  irradiance  Eu(R,zq)  (its  weight)  is  recorded. 
If  the  scalar  irradiance,  in  general  defined  according  to 


E0(  r)  =  L(r,i)dO(i), 
.  11 


(21) 


where  the  integration  is  taken  over  4ir  Sr,  is  desired  at 
the  surface,  i.e.,  £0(R.zo)>  the  photon’s  weight  divided 
by  the  cosine  of  the  angle  its  path  makes  with  the 
vertical  is  recorded,  on  impact  with  the  surface.  If  the 
photon  penetrates  the  surface,  it  also  makes  a  contri¬ 
bution  to  the  associated  irradiance  above  the  surface. 


B.  Validation  Tests  of  the  Monte  Carlo  Code 

The  reciprocity  principle1  can  be  used  to  provide 
some  information  concerning  the  validity  of  the  result¬ 
ing  Monte  Carlo  cbde.  For  example,  let  £%(R,z q)  be 
the  upward  irradiance  exiting  the  sea  surface  due  to  a 
point  source  at  z0  emitting  one  photon  per  second  (Fig. 
2),  and  Eo(zn)  be  the  scalar  irradiance  at  z o  resulting 
from  an  incident  uniform  radiance  distribution  of  unit 
irradiance.  Then  Using  the  reciprocity  principle  it  can 
be  shown  that  (see  Ref.  16  for  the  derivation  of  a 
similar  relationship) 
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E0(z0 )  -  8 rm!  f  El(R,z0)dR.  (22) 

Jo 

This  relationship  was  tested  by  evaluating  the  integral 
using  the  present  code  and  Eo(zo)  using  an  existing  and 
well-tested  Monte  Carlo  code  with  diffuse  illumina¬ 
tion.  The  phase  functions  used  were  similar  to  those 
measured  in  the  ocean,  and  the  computations  agree  to 
within  V2%. 

Unfortunately,  reciprocity  cannot  guarantee  the 
correct  functional  dependence  of  the  point  source  light 
field  on  R  at  the  surface.  This  requires  comparison  of 
the  code’s  output  with  the  results  of  exact  computa¬ 
tions.  For  the  case  of  pure  absorption,  i.e.,  b  =  0,  the 
exact  computation  is  trivial  to  perform,  and  the  Monte 
Carlo  code  agrees  well  with  it.  When  scattering  is 
included,  the  only  point  source  case  for  which  an  ex¬ 
act17  solution  is  available  is  that  of  an  isolated  source  in 
an  infinite  medium,  which  scatters  isotropically,  i.e., 
P(l'  -*•  |)  =  l/4x.  Unfortunately,  the  medium  for 
which  the  present  code  is  applicable  is  semi-infinite 
with  the  detector  at  the  surface.  To  our  knowledge  the 
RTE  has  never  been  solved  exactly  for  such  a  geome¬ 
try.  Elliott,19  however,  has  obtained  an  asymptotic 
solution  for  the  scalar  irradiance  to  the  problem  of  a 
point  in  an  isotropically  scattering  half-space  of  unit 
refractive  index.  Elliott’s  result  for  the  scalar  irradi¬ 
ance  at  the  surface  due  to  a  point  source  emitting  unit 
power  isotropically  at  depth  z0  is 

[E0(rfi,r,)]„  =■  (l  +  e*P<-Vft>.  (23) 

where  tz  -  cz0,  rR  -  cR,  and  k0  and  SM^)  are  constants 
[fe0  =  fco(w0))  which  have  been  tabulated  by  various 
authors.2021  The  asymptotic  solution  is  valid  for  zJR 
«  1.  The  error  is  of  the  order  of  Zq/R5.  The  notation 
[E0(rfi,r2)]iv  means  the  scalar  irradiance  in  nondimen- 
sional  or  scaled  optical  units.  The  actual  scalar  irradi¬ 
ance,  and,  in  general,  the  radiance,  at  the  surface  of  the 
half-space  depends  on  c  as  well  as  zo  and  R,  i.e.,  unlike  a 
homogeneous  ocean  illuminated  by  sunlight,  for  which 
the  radiometric  quantities  depend  on  c  and  depth  only 
through  their  product,  the  point  source  illumination 
produces  a  light  field  that  depends  on  these  quantities 
individually.  The  actual  scalar  irradiance  is  derived 
from  the  scaled  scalar  irradiance  through 


E0{R,z0,c)  =  c^EoItk.t,)]*,  (24) 

where  R  and  z<>  are  111  meters,  and  c  is  in  m_1.  The 
source  power  is  1  W/nm,  and,  since  the  actual  scalar 
irradiance  Ep(R,Zpc)  has  units  of  W/m2nm,  the  scaled 
scalar  irradiance  [2?o(tr,tz]n  has  the  same  units  as  the 
power  (W/nm).  Equation  (24)  is  exact  and  a  manifes¬ 
tation  of  the  inverse  square  decrease  in  the  irradiance 
with  distance  from  the  source.  It  applies  to  other 
radiometric  quantities,  such  as  the  irradiance  and  ra¬ 
diance  as  well. 

The  Monte  Carlo  computations  for  an  isotropically 
scattering  half-space  (with  «o  =  0.9)  illuminated  by  an 
embedded  point  source  are  compared  with  Elliott’s 
asymptotic  theory  in  Fig.  3.  The  individual  panels  in 
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M  (c) 


(b)  (d) 

Fig.  3.  Comparison  between  the  Monte  Carlo  determined  scalar  irradiance  at  the  surface  IE0(tr,ti)]n  in  W/nm,  (noisy  curves)  and  Elliott's 
asymptomatic  solution  for  isotropic  scattering  (smooth  curves),  which  is  valid  for  R  »  zo.  For  all  the  cases  t,  is  0.5,  while  the  physical  depth  of 

the  source  is  varied  from  0.125  m  (a)  to  1  m  (d). 


Fig.  4.  Effect  of  changing  the  refractive  index  of  the  medium  m 
from  1  OOOjFig.  3(c))  to  1.333  for  isotropic  scattering.  Thenoisyline 
is  the  Monte  Carlo  result  for  the  scalar  irradiance  just  beneath  the 
surface  with  m  *  1.333,  and  the  smooth  curve  is  Elliott’s  asymptom¬ 
atic  theory  for  m  =  1.000. 


this  figure  all  refer  to  a  source  at  an  optical  depth  r2  of 
one-half  but  physical  depths  z0  of  1. 0.5, 0.25,  and  0.125 
m.  The  agreement  is  excellent  for  R  »  zo  and  provides 
a  quantitative  measure  of  the  range  of  parameters  over 
which  Elliott's  theory  is  valid,  i.e.,  /?/z0  £  10  or  tr/t*  > 
10.  Statistical  fluctuations  in  the  Monte  Carlo  esti¬ 
mates  are  clearly  evident  for  small  values  of 
[Eo(tr,t?)]/v.  Figure  4  shows  the  effect  of  changing  the 
refractive  index  of  the  scattering  medium  from  1.00  to 
1.33  for  zo  =*  0.5  m  [compare  with  Fig.  3(c)].  For  this 
particular  case  the  scalar  irradiance  at  the  surface  is 
increased  at  all  depths,  and  for  R  a*  10  m  this  increase 
reaches  a  factor  of  4.  This  increase  is  easy  to  under¬ 
stand:  isotropic  scattering  with  large  values  of  «o 
produces  a  very  diffuse  light  field,  for  which  about  half 
of  the  flux  incident  on  the  interface  is  reflected  back 
into  the  medium  and  then  can  be  scattered  back  to  the 
interface. 

The  above  tests  of  the  present  code  suggest  it  is 
correct  and  can  be  used  to  compute  the  irradiance 
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RAYLEIGH  SCATTERING:  cj0  =  0.7;  t,  =  I,  J,  5.  7,  9,  12 


tR 

(°) 

RAYLEIGH  SCATTERING:  u ,  =  0.7;  r,  =■  1.  S.  S,  7.  9,  12 


tR 

(b) 

Fig.  5.  Upwelling  irradiance  just  beneath  the  surface  (£u(fR.r,)].v 
in  W/nm  for  source  optical  depths  1,  3,  5,  7, 9,  and  12,  u>o  of  0-7,  and 
the  Rayleigh  scattering  phase  function:  (alOSr#<  10;  (b)  0  S  tr  < 

30, 


distribution  from  a  point  source  embedded  in  the  wa¬ 
ter. 

C.  Examples  of  the  Computations 
Computations  of  the  distribution  of  upwelling  irra¬ 
diance,  [Eu(TfilT,)]/vor  2?u(fi,z0,c)  (not  the  scalar  irradi¬ 
ance  used  for  comparisons  with  Elliott’s  computa¬ 
tions),  at  the  sea  surface  (just  beneath  the  surface)  has 
been  carried  out  for  a  range  of  values  of  the  parameters 
involved.  Figures  5  and  6  contrast  the  scaled  irradi¬ 
ance  distributions  [Eu(r«,re)]^f  i.e.,  Eu(R^,c)/c2  in  W/ 
nm,  obtained  in  an  ocean  exhibiting,  respectively,  pure 
Rayleigh  scattering  and  scattering  according  to  Pet- 
zold’s  mean  particle  phase  function.  A  point  source 
emitting  a  power  of  1  W/nm  is  placed  at  optical  depths 
of  1,  3,  5,  7,  9,  and  12,  and  o>o  for  the  medium  is  0.7. 
Again,  the  noisy  nature  of  the  results  is  due  to  statisti¬ 
cal  fluctuations  in  the  Monte  Carlo  estimates.  The 
contrast  between  Rayleigh  scattering  (Fig.  5)  and 
strong  forward  scattering  (Fig.  6)  is  seen  to  be  very 
large,  except  for  small  source  depths  and  small  values 


STRONG  rORWARO  SCATTERING:  u,  =  0.7;  Tt  •  I.  J,  5.  7.  9.  12 
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STRONG  rORWARO  SCATTERING:  u,  »  0.7;  T,  »  1 .  J.  5.  7.  9.  1  2 
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Fig.  6.  Upwelling  irradiance  just  beneath  the  surface  [£u(rR,r,)]jv 
in  W/nm  for  source  optical  depths  1, 3,  5, 7. 9,  and  12,  ui0of  0.7,  and  a 
strongly  forward  scattering  phase  function  (a  combination  of  Ray¬ 
leigh  scattering  and  Petzold’s  particle  phase  function  with  cp/cw  » 
7):  (a)  0  <  r«  <  10;  (b)  0  <  r„  <  30. 


of  tr  i.e.,  when  the  source  is  close  to  the  observation 
point.  In  this  regime  the  variation  in  the  irradiance 
with  tr  is  governed  largely  by  geometry  rather  than  by 
the  optical  properties  of  the  water.  Most  photons 
propagate  directly  from  the  source  to  the  detector  with 
no  interaction  with  the  medium.  Scattering  plays  a 
very  small  role  in  determining  the  irradiance.  The 
irradiance  is  proportional  to  z0/X3,  where  X2  =  R2  +  z§ 
(see  below). 

For  larger  values  of  tr  (but  still  small  source  depths) 
the  disparity  between  the  irradiances  for  these  two 
cases  becomes  apparent  with  the  strong  forward  scat¬ 
tering  irradiance  becoming  an  order  of  magnitude 
higher  than  the  Rayleigh  irradiance  by  tr  =  10  for  r,  «* 
1  or  3.  The  rate  of  decay  of  the  irradiance  with  tr 
appears  to  be  independent  of  the  source  depth  for 
sufficiently  large  tr.  This  is  particularly  evident  in 
Fig.  6,  in  which  the  curves  for  various  tr  values  tend  to 
become  parallel  for  large  tr.  Elliott’s  asymptotic  the¬ 
ory  for  isotropic  scattering  predicts  that  this  should  be 
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true  for  the  scalar  irradiance,  since  the  dependence  of 
Eo  on  r.  is  completely  contained  in  the  multiplicative 
'I'o(rJ)  term.  Figure  6  confirms  that  this  result  is  valid 
for  a  strongly  forward  scattering  phase  function  as 
well.  This  is  the  diffusion  regime;  the  light  field  is 
completely  dominated  by  multiple  scattering.  Pho¬ 
tons  scatter  many  times  before  encountering  the  sur¬ 
face.  In  the  diffusion  approximation  the  irradiance 
would  be  expected  to  vary  with  R  according  to 

EJK,zt„c)  i 

- exp(-*dTV/c),  (25) 

f  rx 

where 

Tx  +  rl- 

and  kj  is  explicitly  provided  in  terms  of  a>0  and  the 
phase  function  asymmetry  parameter  (or  average  co¬ 
sine)  g  through-1 

kd  , _ _ . _ 

y  =  -w„)(l  -uy>).  (26) 

In  this  approximation,  kd  is  the  downwelling  (and  up- 
welling)  irradiance  attenuation  coefficient  ( Kd )  of  the 
light  field  when  illuminated  by  solar  radiation.  The 
diffusion  approximation  is  valid  only  at  very  large 
values  of  rR  and  rz  and  values  of  w0  near  l,i.e.,  hjc  «  1. 
The  irradiance  reaching  the  surface  is  significantly 
larger  with  strong  forward  scattering  than  with  Ray¬ 
leigh  scattering,  and  this  phase  function  effect,  ex¬ 
pressed  through  the  parameter  g  in  Eq.  (26),  becomes 
increasingly  stronger  as  tr  and  rz  increase.  For  large 
tr  the  exponential  dominates  the  algebraic  function, 
and  the  result  is  essentially  exponential  decay.  It  is 
natural  to  expect  the  exponential  decay  constant  * 
associated  with  this  exponential  decay  to  be  related  to 
the  attenuation  coefficient  of  irradiance  K»  in  the 
asymptotic  light  regime.22  The  asymptotic  decay  co¬ 
efficients  Km  have  been  computed  for  the  two  cases 
examined,  along  with  *,  which  was  determined  by  fit¬ 
ting  the  computation  for  large  tp  (tr  >  10)  to  Eq.  (25) 
with  kd  replaced  by  k.  These  are  compared  in  Table 
IV,  which  gives  the  results  for  Rayleigh  and  forward 
scattering  with  wo  a  0.7  (first  two  rows)  and  for  forward 
scattering  with  wo  =  0.9  (last  row).  Figure  5  shows  that 
the  Rayleigh  irradi4nce  is  very  noisy  for  large  tz  and 
large  tr,  and  hence  the  k/c  determined  from  the  com- 
puta.ions  is  likely  to  be  inaccurate.  The  strong  for¬ 
ward  scattering  irradiance,  on  the  other  hand,  (Fig.  6), 
is  much  better  behaved  at  large  tr,  and  we  can  expect  a 
more  accurate  k/c.  The  values  of  k/c  for  wo  =  0.7  (first 
two  rows  in  Table  IV)  usually  agree  with  K„  within  the 
accuracy  of  the  k/c  determination,  with  the  exception 
of  the  Rayleigh  scattering  case  with  tz  =  1  and  3.  The 
case  of  strong  forward  scattering  with  w0  =  0.9  (last  row 
in  Table  IV)  was  included  for  two  reasons:  the  irradi¬ 
ance  decays  much  more  slowly  with  rR,  enabling  larger 
values  of  tr  to  be  reached  in  the  computations  (the  fits 
for  k/c  include  only  results  for  tr  >  15  rather  than  the 
•  r  >  10  in  the  other  chses);  and  the  diffusion  approxi¬ 
mation  becomes  more  appropriate  for  larger  values  of 
w().  The  agreement  between  k/c  and  is  seen  to  be 


Table  IV.  Darloea  Valuaa  o I  k/c  lot  Comparison  with  K./c 


P(0) 

wo 

czn 

KJc 

1 

3 

5 

7 

9 

12 

Rayleigh 

0.7 

0.929 

1.050 

0.864 

0.939 

1.013 

1.010 

0.821 

Forward 

0.7 

0.470 

0.423 

0.449 

0.465 

0.480 

0.414 

0.467 

Forward 

0.9 

0.245 

0.226 

0.232 

0.238 

0.242 

0.240 

0.253 

PIGMENT  CONCENTRATION  0.10  mj/m'  I,  t  10,  JO.  50,  70,  90m 


T» 

(o) 

PIGMENT  CONCENTRATION  0.50  2.  =  10.  30.  50.  70.  90m 


tR 

(b) 

Fig.  7.  Vpwelling  irradiance  just  beneath  the  surface 
in  W/nm  and  source  physical  depths  10,  30,  50,  70,  and  90  m  for 
pigment  concentrations  of  0.1  (a)  and  0.5  (b)  mg/m3.  Note:  to 
compute  the  actual  irradiance  (W/m!nro)  the  scaled  irradiance  given 
here  must  be  multiplied  by  c2. 

better  and  more  stable  with  source  depth  for  w0  «*  0  9 
compared  with  the  other  case.  In  the  two  forward 
scattering  cases,  which  provide  a  realistic  approxima¬ 
tion  to  the  real  ocean,  equality  between  k/c  and  is 
obtained  within  10%  for  all  source  depths. 

Figure  7  shows  the  results  of  computations  of 
lEu(rR,r*)]w  carried  out  for  source  depths  of  10, 30,  50, 
7H  rH  90  m  in  an  ocean  exhibiting  optical  properties 
at  480  nm  consistent  with  the  bio-optical  model  de¬ 
scribed  above  for  pigment  concentrations  of  0. 1  and  0.5 
mg/m3.  Such  concentrations  are  typical  of  central 
ocean  gyres  at  times  of  medium  to  high  productivity. 
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These  results  are  very  similar  to  the  forward  scattering 
cases  described  above.  As  expected,  the  higher  pig¬ 
ment  concentration  results  in  a  lower  level  of  irradi- 
ance  at  the  surface  for  all  source  depths.  For  a  de¬ 
tailed  comparison  of  these  two  cases,  one  must  work 
with  actual  (as  opposed  to  scaled)  irradiances.  As 
described  above,  the  actual  irradiance  is  derived  from 
the  scaled  irradiance  by  multiplication  by  c2.  For  C  of 
0.1  and  0.5  mg/m1,  c  is  0.12  and  0.28  m_I,  respectively. 
Thus  to  compare  the  two  cases  at,  say  R  =  30  m,  we 
must  compute  j-r  =  cR  and  multiply  the  nondimen- 
sional  irradiance  by  c2.  This  gives  at  z0  =  10  m  irra¬ 
diances  of  ~1.15  X  10~5  and  6.3  X  10~6  W/m2  nm  for  C 
of  0.1  and  0.5  mg/m1,  respectively,  if  the  source  power 
is  1  W/jim.  It  may  seem  surprising  at  first  that  the 
surface  irradiance  in  the  clearer  water  is  only  1.8  times 
greater  than  that  in  the  more  turbid  water.  However, 
one  must  remember  that  the  plankton  scatter  much 
more  strongly  than  they  absorb  (ojq  =  0.88)  and  that 
most  of  the  scattering  is  in  the  naar  forward  direction 
(Fig.  1 ).  Hence,  in  the  first  approximation,  the  irradi¬ 
ance  decay  is  governed  by  the  absorption  and  geome¬ 
try.  At  480  nm  the  absorption  coefficients  for  the 
plankton  and  detritus  are  0.012  and  0.032  m-1,  respec¬ 
tively,  for  the  low  and  high  pigment  concentrations. 
F rom  the  geometry  of  the  problem  the  irradiance  must 
travel  at  least  31.6  m  before  reaching  the  surface;  thus 
considering  only  the  differential  absorption  loss,  the 
ratio  of  the  two  irradiances  should  be  exp[-31  6(0.012 
-  0.032)]  =  1.88,  in  good  agreement  with  the  result 
derived  from  Fig.  7. 

D.  Analytic  Approximations  to  the  Monte  Carlo  Results 

It  would  be  useful  to  be  able  to  describe  the  irradi¬ 
ance  distribution  on  the  surface  with  simple  analytic 
formulas.  The  diffusion  approximation  above,  and 
Elliott’s  result,  are  examples  for  which  this  is  possible, 
their  value  is  limited  since  they  basically  apply  in  a 
regime  (rR  »  1)  where  the  irradiance  is  very  small  or 
negligible.  What  is  really  needed  is  to  be  able  to 
describe  [£u(r/?,r,)J  v  for  the  more  important  region  0 
<  th  <  10,  where  most  of  the  irradiance  is  found.  The 
success  of  the  simple  computation  above  in  explaining 
the  Monte  Carlo  result  for  small  zq  suggests  a  simple 
model  to  explain  quantitatively  the  results  in  this  in¬ 
terval.  Consider  the  irradiance  on  the  surface  due  to 
photons  emitted  from  a  point  source  of  unit  power  at  a 
depth  2d  in  the  absence  of  scattering.  This  is 

F.JR,zn.c)  =  cos  '#  exp<  -azg/cosfi),  (27) 

where  tan#  =  R/zn.  Rewriting, 

K,(  R.z„.c)  i  t. 

,  **  |£u(rR.r,)JiV  «  ■—  -j  exp(-arv/c),  (28) 

c-  4*  r'x 

where  a  =  c(r  y  was  defined  earlier).  Now  in  the  case 
of  a  strongly  forward  scattering  medium,  most  photons 
are  scattered  with  very  little  change  in  direction,  so  we 
expect  this  to  still  be  approximately  correct,  but  with  a 
replaced  by  an  effective  attenuation  coefficient  k, 
which  is  expected  to  be  related  to  Kd  in  some  manner. 


Note,  however,  the  Kd  is  not  a  uniquely  defined  con¬ 
stant;  in  general,  it  is  a  function  of  depth. 

There  are  probably  many  ways  to  determine  a  suit¬ 
able  value  of  k,  for  example,  one  could  fit  the  Monte 
Carlo  results  to 


£u<R.*o.c>  .  . 

- -J - = 


y  '  exp(— Arx/c) 
4x 


(28') 


and  determine  A  and  k.  Although  this  technique  pro¬ 
vides  excellent  fits  to  the  Monte  Carlo  results  in  the 
nondiffusion  regime,  it  is  not  easy  to  interpret  A  and  k 
in  terms  of  physically  measurable  quantities,  e  g.,  c, 
Kd-  The  approach  taken  here  is  based  on  trying  to  fit 
the  total  power  reaching  the  surface  to  Eq.  (28)  with  a 
replaced  by  k.  The  total  power  reaching  the  surface  is 
found  by  multiplying  Eq.  (28)  by  2nc‘1RdR  and  inte¬ 
grating  from  0  to  ®,  i.e., 


P10U i  =  f  2 *EJR,z0,c)RdR.  (29) 

JO 

This  integral  can  be  evaluated  in  terms  of  an  exponen¬ 
tial  integral  yielding  the  remarkably  simple  formula 


P.uiai  =  \  E2(~  =  2  E2{hz0),  (30) 

where 


which  is  tabulated  in  Abramowitz  and  Stegun.22  The 
total  power  reaching  the  surface  thus  depends  only  on 
k  and  the  depth  of  the  source  z0. 

The  determination  of  k/c  was  made  by  trial  and 
error  with  the  restriction  that  only  computations  for 
which  PtoUi  >  10-5  would  be  included  for  a  point  source 
of  unit  power.  This  restriction  was  placed  on  the 
analysis  because  (1)  it  encompasses  the  most  interest¬ 
ing  region  and  (2)  it  was  believed  that  smaller  values  of 
Ptoui  would  place  the  problem  in  the  diffusion  regime 
for  which  even  the  form  of  Eq.  (28)  is  incorrect.  Brief¬ 
ly,  in  the  first  attempt,  the  traditional  value  of  Kd  at 
the  surface,  i.e.,  a  +  bt,,  was  tried  for  k.  This  value 
worked  well  for  larger  values  of  PtoUi,  i.e.,  smaller  val¬ 
ues  of  r2,  but  overestimated  PtoU|  for  larger  values  of  r2. 
Realizing  that  Kd  varies  with  depth,  the  next  try  was  to 
estimate  k  by  the  asymptotic  value  of  Kd,  which  was 
computed  for  each  simulation.  This  time  good  agree¬ 
ment  between  the  Monte  Carlo  results  and  Eq.  (30) 
was  found  at  large  values  and  at  very  small  values  of 
ftoui.  but  for  intermediate  values,  Ptotai  *  10~s,  it  was 
underestimated.  A  linear  combination  of  these  two  k 
values  did  not  yield  better  results.  Finally,  it  was 
decided  to  try  a  linear  combination  weighted  by  a 
function  of  r,  so  that  for  small  values  of  r2  the  value  of 
Kd  at  the  surface,  iC8U rf,  would  be  used,  while  for  very 
large  values  of  t2  the  asymptotic  Kd  value  A",  would  be 
used.  This  scheme  worked  fairly  well  for  10~f'  <  Pt„,a| 
<  1 ;  however,  it  was  felt  to  be  unsatisfactory  because, 
although  Ksurf  could  be  easily  measured  in  the  ocean, 
measurement  of  K .  is  difficult.  This  blemish  is  easy 
to  remove,  because  K9l,rf  and  K„  are  found  to  be  closely 
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Mon*e  Carlo  Total  Powar  (Walls/nm) 

Fig.  8.  Comparison  between  the  total  power  just  beneath  the  sea 
surface  (from  a  source  of  unit  power)  computed  via  Monte  Carlo 
techniques  with  that  computed  using  the  analytic  model  with  k 
determined  from  Eqs.  (31)-(33). 


related  in  the  present  model  for  the  pigment  range  0.01 
<  C  <  4.4  mg/m3: 

K.  K  , 

0.0836  < - -  <  0.114. 

c  c 

Similar  relationships  have  been  suggested  by  other 
investigators.24  25  By  trial  and  error,  the  final  rela¬ 
tionship  for  K „  was  taken  to  be 

«-  =  Km„  +  0.09c,  (31) 

and  the  best  value  of  k  to  be  used  in  Eq.  (30)  is 

fe-7(r«)K.urf+ll-7MK.,  (32) 

where  the  weighting  function  7  (r2)  is  given  by 

t(t,)  =  exp^-0.08-K^  =  exp(-0.08K„urfz0).  (33) 

Figure  8  compares  the  Monte  Carlo  computations  of 
the  total  power  with  that  computed  from  Eq.  (30) 
(Analytic  Total  Power)  with  k  determined  from  Eqs. 

(31 ) — (33)  for  pigment  concentrations  ranging  from 

0.01  to  4.4  mg/m3  and  source  depths  from  10  to  90  m. 
Ksur f  is  taken  from  Eq.  (16c).  Equation  (30)  clearly 
provides  an  excellent  approximation  to  the  total  power 
reaching  the  surface  from  a  point  source  of  unit  power 
for  10~5  <  P total  :£  1.  and  Ks urf  can  be  directly 

related  to  C  using  Eqs.  (13)  and  (16),  i.e.,  Ksur{  =  0.020 
+  0.052C06,  and  =  0.022  +  0.087 C06.  These  pro¬ 
vide  k  as  a  function  of  C  and  z0  directly  through  Eqs. 

(32)  and  (33).  . 

Having  determined  the  value  of  k  that  reproduces 
the  total  power,  this  value  can  be  substituted  into  Eq. 
(28)  (replacing  ri)  to  estimate  the  distribution  of  irradi- 
ance  on  the  sea  surface  [Eu(th,t2)]^.  These  distribu¬ 
tions  for  tr  <  10  are  shown  in  Fig.  9  with  pigment 
concentrations  of  0.01, 0.5, 0.1,  and  0.5  mg/m3.  Clear¬ 
ly,  the  simple  model  reproduces  the  Monte  Carlo  com¬ 
putations  in  this  nondiffusion  regime  very  well.  It  is 
interesting  to  note  that  these  results  imply  that  an 
estimate  of  Kj  can  be  made  at  night  by  releasing  a 


point  source  of  light  in  the  water,  measuring  the  irradi- 
ance  at  the  surface  as  it  sinks,  i.e.,  as  a  function  of  the 
depth  of  the  source,  and  fitting  the  results  to  the  equa¬ 
tions  developed  here  to  determine  k. 

This  simple  model  can  be  used  to  estimate  the  sur¬ 
face  manifestation  of  bioluminescent  emissions,  given 
the  optical  properties  of  the  medium,  or  as  seen  in  the 
present  model,  the  pigment  concentration.  For  exam¬ 
ple,  one  can  determine  the  fraction  of  the  total  power 
reaching  the  surface  in  a  circle  of  radius  R,  i.e.,  the 
effective  size  of  the  spot  on  the  sea  surface.  The  power 
reaching  the  surface  enclosed  in  a  circle  of  radius  R  is 
given  by 

P(R)  =  f  2*EJR',zu)R'dR'.  (34) 

Jo 

This  is  easily  evaluated  in  terms  of  exponential  inte¬ 
grals,  and  the  result  is 

P(R)  *  V,{£:2(Az0)  -  -  £2(Ar0n)].  (35) 

V 

where  t?  =  vl  +  P2/z2.  The  fraction  / of  the  total  power 
contained  within  R  is  /  =  P(R/)/Plo ui,  where  Rf  is  the 
value  of  R  corresponding  to  f.  Thus  the  radius  of  the 
circle  containing  the  fraction  /of  the  total  power  reach¬ 
ing  the  surface  is  found  by  solving  the  equation. 

Et{kt0r,)  =  „(1  -  f)E2(fez0)  (36) 

for  ij.  The  resulting  values  of  Rf  for  /  =  0.5  and  0.9  are 
presented  in  Fig.  10.  Choosing  /  =  0.9  gives  the  effec¬ 
tive  size  of  the  spot  on  the  sea  surface.  The  accuracy  of 
this  determination  can  be  estimated  from  Fig.  11  in 
which  TRf  (=  cR/)  determined  from  the  individual  Mon¬ 
te  Carlo  simulations  is  compared  with  the  same  quan¬ 
tity  determined  analytically  using  Eqs.  (31^— (36). 
The  excellent  agreement  for  {  =  0.5  is  seen  to  be  de¬ 
graded  only  slightly  for  /  =  0.9.  The  possibility  of 
inverting  Eq.  (36)  for  z0  given  R  and  some  measure  of  k 
is  discussed  in  the  Appendix. 

Since  our  computations  of  Eu(R^,c)  are  for  a  point 
source,  the  results  represent  the  Green’s  function  for 
the  time-independent  transport  equation  and,  there¬ 
fore,  Eq.  (28),  with  a  replaced  by  k,  represents  an 
approximation  to  the  Green’s  function.  Thus,  for  an 
extended  distribution  of  sources  of  power  P(ro)  dis¬ 
tributed  throughout  the  medium,  the  upwelling  irradi- 
ance  just  beneath  the  surface  will  be  given  by 


X  exp[— *|  r0  -  r,|  )dr0,  (37) 

where  rs  is  the  point  on  the  surface  at  which  the  irradi- 
ance  is  required,  and  d ro  is  an  abbreviation  for 
dxodyodzo. 

Finally,  the  irradiance  just  above  the  sea  surface,  i.e., 
the  component  of  (Eu(T/e,r2)]^  that  is  transmitted 
through  the  ai.*-sea  interface  remains  to  be  discussed. 
This  quantity  is  easy  to  derive  by  combining  the  radi¬ 
ance  just  beneath  the  sea  surface  with  the  Fresnel 
equations  for  specular  reflection  and  is  carried  out  in 
the  Monte  Carlo  code;  however,  these  computations 
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PIGMEN  T  CONCENTRATION  0.01  mg/m5;  i,  =  10,  30.  50.  70.  90m  PIGMENT  CONCENTRATION  0.10  mg/m*;  «.  ■  10.  30,  50.  70.  90m 


Fig.  9.  Upwelling  irradiance  just  beneath  the  surface,  [£u(rR,r;)]M  in  W/nm,  and  source  physical  depths  10,  30,  50,  70,  and  90  m  for  pigment 
concentrations  of  0.01  (a),  0.05  (b),  0.10(c),  and  0.50  (d)  mg/m3.  The  noisy  curves  are  the  Monte  Carlo  results,  and  the  smooth  curves  are  com¬ 
puted  from  the  analytical  model. 


Fig  10.  Analytic  approximation  to  the  radius  of  the  circle  at  the  sea 
surface  /?/  containing  a  fraction  /  of  the  total  power  plotted  as  a 
function  of  the  source  depth 


are  less  important  than  those  for  [£u(tr»’-2)]iv»  since  the 
actual  irradiance  leaving  the  surface  will  depend  sig¬ 
nificantly  on  the  roughness  of  the  surface,  which  here 
has  been  assumed  to  be  perfectly  flat.  The  results  are 
similar  to  [£u(tr,Tj)];v1  however,  when  r,  is  small  the 
irradiance  shows  a  rapid  decrease  with  tr,  beyond  tr  « 
r„  because  unscattered  photons,  and  photons  scat¬ 
tered  through  very  small  angles,  strike  the  sea  surface 
at  incident  angles  greater  than  the  critical  angle  for 
total  internal  reflection.  It  is  of  interest  to  see  how 
well  the  analytical  model  above  represents  the  irradi¬ 
ance  above  the  surface.  This  is  shown  in  Fig.  12  (com¬ 
pare  with  Fig.  9  for  [£u(rR,Tj)^|,  in  which  the  scaled 
irradiance  is  plotted  as  a  function  of  tr  and  Zo  for 
pigment  concentrations  ranging  from  0.01  to  0.5  mg/ 
m3  along  with  the  predictions  of  the  analytical  model 
derived  by  multiplying  Eq.  (28),  with  a  replaced  by  k, 
by  the  Fresnel  transmittance  evaluated  for  an  incident 
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Fig.  11.  Comparison  between  the  analytic  approximation  to  the  radius  of  the  circle  at  the  sea  surface,  tR/  e  cR/t  containing  a  fraction  /  of  the 
total  power,  plotted  as  a  function  of  the  source  optical  depth  r,  a  c2q :  (a)  /  =  50%;  (b)  /  =  90%. 


PICUENT  CONCENTRATION  0.01  mg/m  ;  i  «  10.  30.  SO,  70,  90m 


PIGMENT  CONCENTRATION  0.10  mg/m  :  I.  =  10,  30.  SO.  70.  90m 


PIGMENT  CONCENTRATION  0.05  mg/m';  =  10.  30,  50.  70,  90m  1 


PIGMENT  CONCENTRATION  0.50  mg/m  ;  =  10.  30.  50,  70.  90m 


Fig.  12  Scaled  upwelling  irradiance  W/nm  just  above  the  surface  for  source  physical  depths  10,  30,  50,  70,  and  90  m  and  pigment 
concentrations  of  0.01  (a),  0.05  (b),  0. 10  (c),  and  0.50  (d)  mg/m3.  The  noisy  curves  are  the  Monte  Carlo  results,  and  the  smooth  curves  are  com¬ 
puted  from  the  analytical  model  accounting  for  refraction  and  reflection  at  the  sea-air  interface. 
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angle#,  where  tan0  =  tr/t:.  The  rapid  reduction  of  the 
irradiance  with  tr,  after  tr  ~  tz,  is  clearly  evident  for 
the  smaller  values  of  z0.  The  analytical  model,  of 
course,  predicts  that  the  irradiance  above  the  surface  is 
identically  zero  for  incident  angles  greater  than  6  * 
48.6°,  i.e.,  the  vertical  lines  in  Fig.  12.  Clearly  the 
analytical  model  fi.ts  the  Monte  Carlo  data  well  for 
values  of  tr  and  tz  such  that  6  <  48.6°,  and  for  small 
values  of  zq,  this  regime  contains  most  of  the  total 
power  exiting  the  surface.  For  large  values  of  z0  the 
light  field  at  the  surface  becomes  very  diffuse  and  the 
effect  of  the  critical  angle  is  hardly  observed.  In  this 
regime,  a  significant  fraction  of  the  total  power  exiting 
the  water  is  outside  the  range  of  applicability  of  the 
analytic  model.  Figure  13  shows  the  comparison  be¬ 
tween  the  total  power  computed  using  the  analytical 
model  and  the  Monte  Carlo  computations.  Clearly 
the  analytical  model  provides  excellent  results  for 
cases  where  the  total  power  exiting  the  surface  is  great¬ 
er  than  ~10~2  of  the  source  power,  but  for  the  10'5- 
10~3  range  it  underestimates  the  total  power  by  as 
much  as  a  factor  of  2. 

IV.  Summary  and  Conclusions 

A  model  of  the  optical  properties  of  the  ocean,  pro¬ 
viding  the  absorption  and  scattering  coefficients  of  the 
medium  as  nonlinear  functions  of  the  concentration  of 
pigments  associated  with  phytoplankton  and  their  im¬ 
mediate  detrital  material,  is  presented.  Monte  Carlo 
computations  of  the  attenuation  coefficient  of  down- 
welling  irradiance  Kd  for  an  ocean-atmosphere  system 
illuminated  by  the  sun  at  zenith  agree  well  with  experi¬ 
mental  data  and  demonstrate  the  validity  of  such  a 
model  for  studying  the  influence  of  biomass  on  the 
propagation  to  the  surface  of  light  generated  through 
bioluminescence.  The  radiative  transfer  equation  for 
the  irradiance  at  the  sea  surface  resulting  from  illumi¬ 
nation  by  a  point  source  embedded  in  the  water  is 
solved  by  Monte  Carlo  techniques.  The  solution  tech¬ 
nique  is  validated  through  comparison  with  an  asymp¬ 
totic  analytic  solution  for  isotropic  scattering.  The 
computations  show  that  the  irradiance  distribution 
just  beneath  the  surface  as  a  function  of  R,  the  distance 
measured  along  the  surface  from  a  point  on  the  surface 
vertically  above  the  source,  is  described  by  two  re¬ 
gimes:  ( 1 )  a  regime  in  which  the  irradiance  is  governed 
mostly  by  absorption  and  geometry,  with  scattering 
playing  a  negligible  role,  the  near  field;  and  (2)  a  regime 
in  which  the  light  field  at  the  surface  is  very  diffuse  and 
the  irradiance  decays  approximately  exponentially  in. 
R,  with  a  decay  coefficient  and  is  a  very  weak 
function  of  the  source  depth,  the  diffusion  regime. 
The  near  field  is  of  primary  interest  because  it  contains 
most  of  the  power  reaching  the  sea  surface. 

An  analytical  model  of  the  irradiance  distribution 
just  beneath  the  surface  as  a  function  of  R,  the  source 
depth,  and  the  pigment  concentration  for  the  near 
field  is  presented.  This  model  is  based  on  the  observa¬ 
tion  that  at  most  scattering  events  the  change  in  the 
photon’s  direction  is  slight  and,  therefore,  ineffective 
in  attenuating  the  irradiance.  A  solution  for  the  irra- 


Monte  Carlo  Total  Power  (Watt»/nm) 

Fig.  13.  Comparison  between  the  total  power  just  above  the  sea 
surface  (from  a  source  of  unit  power)  computed  via  Monte  Carlo 
techniques  with  that  computed  using  the  analytic  model  with  k 
determined  from  Eqs.  (31  )-(33). 

diance  from  the  point  source,  then,  is  first  carried  out 
ignoring  scattering  altogether;  however,  recognizing 
that  backscattering  will  attenuate  the  irradiance,  the 
absorption  coefficient  is  replaced  by  an  effective  atten¬ 
uation  coefficient  k.  This  effective  attenuation  coeffi¬ 
cient  is  determined  by  fitting  the  total  power  just 
beneath  the  surface  determined  from  the  Monte  Carlo 
computations  to  the  analytical  model.  The  resulting  k 
is  closely  related  to  Kdt  with  k  «  a  +  bb  for  small  source 
depths,  and  k  K„  as  the  source  depth  becomes  very 
large.  A  relationship  is  developed  giving  k  as  a  func¬ 
tion  of  the  source  depth  and  the  pigment  concentra¬ 
tion,  which  reproduces  the  total  power  incident  on  the 
surface  with  remarkable  precision.  Using  k  deter¬ 
mined  in  this  manner,  the  Monte  Carlo  irradiance  as  a 
function  of  R  and  source  depth  in  the  near-field  regime 
can  be  approximated  with  high  accuracy.  These  re¬ 
sults  indicate  that  Kd  can  be  estimated  at  night  by 
releasing  a  point  source  in  the  water,  measuring  the 
irradiance  at  the  surface  as  it  sinks,  and  fitting  the 
measurements  to  the  relationships  developed  here  to 
determine  k .  The  analytic  model  also  enables  estima¬ 
tion  of  the  source  depth  and  power  from  the  irradiance 
distribution  just  beneath  the  surface. 

The  irradiance  exiting  the  water  was  also  deter¬ 
mined  in  the  Monte  Carlo  simulations,  but  in  this  case 
the  usefulness  of  the  analytical  model  is  limited  be¬ 
cause  it  predicts  that  the  exiting  irradiance  is  zero  for 
values  of  R  greater  than  Rc  =  z0  tan#c,  where  6C  is  the 
critical  angle  for  total  internal  reflection.  When  the 
depth  of  the  source  is  small,  only  an  insignificant 
amount  of  power  falls  outside  a  circle  of  radius  Rc,  but 
when  it  is  large  this  is  not  always  the  case  and  the 
analytic  model  can  lead  to  significant  errors  in  the 
estimation  of  the  total  power  exiting  the  ocean.  How¬ 
ever,  for  R  <  Rc  the  analytical  model  fits  the  Monte 
Carlo  computed  irradiance  above  the  surface  with  ex¬ 
cellent  accuracy. 

This  work  received  support  from  the  Office  of  Naval 
Research  under  contract  N00014-84-K-0451  as  part  of 
the  Biowatt  Program. 
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Mont#  Carlo  tR) 

(b) 

Fig.  14.  Comparison  between  the  analytic  approximation  to  the 
radius  of  the  circle  at  the  sea  surface,  £  cRf,  containing  a  fraction  f 

of  the  total  power,  plotted  as  a  function  of  the  source  optical  depth  r, 
Hcz0.  In  this  case  k  has  been  approximated  by  K,u,f.'  (a)  f  =  50%;  (b) 
/  =  90%. 

Appendix:  Surface  Estimation  of  the  Source  Depth  and 
Power 

Figure  10  suggests  the  possibility  of  remote  estima¬ 
tion  of  the  emission  depth  given  P(R)/Pt0 ui,  i.e.,  if  k 
were  known,  z0,  the  emission  depth,  could  be  deter¬ 
mined  remotely  from  observations  of  the  spot  size  Rf. 
The  difficulty,  however,  is  in  the  determination  of  k, 
since  this  requires  knowing  zo,  i.e.,  Eq.  (33).  Thus  it  is 
of  interest  to  know'  how  well  the  analytic  formulation 
predicts  without  using  z0  to  compute  k.  A  reason¬ 
able  approximation  to  k  is  K„ur f,  which,  although  it 
does  not  produce  as  good  a  fit  to  the  total  power  as  that 
determined  using  ,Eqs.  (31)-(33)  (Fig.  8),  yields  a 
[Eu(tr,t2)]n  as  a  function  of  rfl  and  r2  having  the  cor¬ 
rect  shape.  Figure  14  compares  the  Monte  Carlo  and 
analytic  determinations  of  tr(  using  Eq.  (36)  with  k  = 
K,„rf.  The  results  show  only  a  slight  degradation  from 
those  in  Fig.  11  obtained  using  the  optimum  value  of  k. 
In  both  cases  if  =  50  and  90%)  most  of  the  degradation 
occurs  for  larger  values  of  trs,  for  which  the  approxi¬ 
mate  formulation  ( k  =  ffsurf)  overestimates  the  value  of 
tr,  given  r2.  Figure  15  shows  the  r2  retrievals  from  the 


(o) 


(b) 

Fig.  15.  Comparison  between  the  source  optical  depth,  t,  b  cz0, 
retrieved  by  inverting  Eq.  (36)  for  various  pigment  concentrations 
and  source  physical  depths:  (a)  f  =  50%;  (b)  f  =  90%. 


Monte  Carlo  determined  tr/  derived  by  inverting  Eq. 
(36)  with  ft  =  K% urf.  Clearly,  such  an  inversion  appears 
possible  for  moderate  values  of  r2,  i.e.,  t2  <  10,  even  in 
the  presence  of  the  noise  inherent  in  the  Monte  Carlo 
simulations.  This  inversion  in  the  presence  of  Monte 
Carlo  noise  suggests  that  it  may  be  feasible  in  the  real 
ocean. 

Another  possibility  for  inverting  the  irradiance  dis¬ 
tribution  is  to  use  the  pigment  concentration  to  deter¬ 
mine  If 8Urf  using  Eq.  (16d).  This  would  approximate 
the  Kd  value,  which  would  actually  be  measured  (dur¬ 
ing  the  day  with  the  sun  near  the  zenith)  near  the 
surface.  Figure  16  shows  the  result  of  using  this  Xaurf 
to  retrieve  z<>.  The  retrievals  become  progressively 
poorer  at  all  source  depths  as  the  pigment  concentra¬ 
tion  increases,  i.e.,  the  retrieval  becomes  poorer  as  r2 
increases.  Nevertheless,  the  retrievals  are  remarkably 
accurate  with  -0.15  <  6z0/z0  <  0.08  for  /  =  50%  and 
0.05  £  C  £  1  mg/m3.  This  is  particularly  impressive 
for  the  case  where  z0  =  120  m  and  C  =  1  mg/m3  for 
which  c  e  0.46  m- 1 ,  t2  =  55.2,  and  the  probability  that  a 
photon  will  reach  the  surface  without  interacting  with 
the  medium  is  <10-24!  Note  that  the  value  of  the 
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earn  attenuation  coefficient  c  is  not  used  anywhere  in 
lese  retrievals. 

Once  zq  is  determined,  the  power  emitted  by  the 
Durce  (W/nm)  can  be  estimated  by  adapting  Eq.  (30) 
3  a  source  of  power  Pq,  i.e., 


It  is  unlikely  that  this  inversion  can  be  extended  to 
leasurements  of  the  upwelling  irradiance  above  the 
urface,  because  the  analytic  model  does  not  provide  a 
ood  approximation  to  the  irradiance  at  large  R,  and 
he  distribution  above  the  sea  surface  will  be  a  function 
f  the  surface  roughness. 

tote  added  in  proof:  The  angle  0  shown  in  Fig.  2  is  incorrect.  9  is 
le  angle  between  nadir  and  the  photon  direction,  i.e.,  the  correct 
ngle  0  is  the  supplement  of  the  angle  shown  in  the  figure. 


True  z„  (Maters) 


teferences 

1.  K.  M.  Case,  “Transfer  Problems  and  the  Reciprocity  Principle,” 
Rev.  Mod.  Phys.  29, 651  (1957). 

2.  A.  Morel  and  L.  Prieur,  “Analysis  of  Variations  in  Ocean  Color,” 
Limnol.  Oceanogr.  22,  709  (1977). 

3.  By  the  term  chlorophyll  a  we  mean  the  concentration  (mg/m3)  of 
chlorophyll  a  and  ail  chlorophyll-like  pigments  which  absorb  in 
the  same  spectral  banes  as  chlorophyll  a,  such  as  phaeophytin  a, 
and  are  contained  in  phytoplankton  or  in  their  detrital  materi¬ 
als.  The  sum  of  the  concentrations  of  chlorophyll  a  and  phaeo¬ 
phytin  a  is  frequently  used  as  an  indicator  of  plankton  biomass. 
It  is  usually  referred  to  as  the  pigment  concentration. 

4.  K.  S.  Baker  and  R.  C.  Smith,  “Bio-optical  Classification  and 
Model  of  Natural  Waters.  2,”  Limnol.  Oceanogr.  27, 500  (1982). 

5.  A.  Morel,  “Optical  Properties  of  Pure  Water  and  Pure  Sea 
Water,”  in  Optical  Aspects  of  Oceanography,  N.  G.  Jerlov  and 
E.  Steemann  Nielsen,  Eds.  (Academic,  New  York,  1974). 

6.  A.  Morel,  “In-water  and  Remote  Measurement  of  Ocean  Color,” 
Boundary-Layer  Meteorol.  18, 177  (1980). 

7.  H.  R.  Gordon  and  A.  Y.  Morel,  Remote  Assessment  of  Ocean 
Color  for  Interpretation  of  Visible  Satellite  Imagery:  A  Re¬ 
view  (Springer-Verlag,  New  York,  1983). 

8.  L.  Prieur  and  S.  Sathyendranath,  “An  Optical  Classification  of 
Coastal  and  Oceanic  Waters  Based  on  the  Specific  Absorption  of 
Phytoplankton  Pigments,  Dissolved  Organic  Matter,  and  Other 
Particulate  Materials,”  Limnol.  Oceanogr.  26, 671  (1981). 

9.  R.  C.  Smith  and  K.  S.  Baker,  “The  Bio-optical  State  of  Ocean 
Waters  and  Remote  Sensing,”  Limnol.  Oceanogr.  23, 247  (1978). 

0.  R.  C.  Smith  and  K.  S.  Baker,  “Optical  Classification  of  Natural 
Waters,”  Limnol.  Oceanogr.  23, 260  (1978). 

1.  L.  A.  Hobson,  D.  W.  Menzel,  and  R.  T.  Barber,  “Primary  Pro¬ 
ductivity  and  the  S^zes  of  Pools  of  Organic  Carbon  in  the  Mixed 
Layer  of  the  Ocean,”  Mar.  Biol.  19, 298  (1973). 

2.  S.  Sathyendranath,  “Influence  des  Substances  en  Solution  et  en 
Suspension  dans  lee  Eaux  de  Mer  sur  L’absorption  er  La  Reflec¬ 
tance.  Modelisatian  et  Applications  a’la  Teledetection,”  Ph.D. 
Thesis,  3rd  cycle,  U.  Pierre  et  Marie  Curie,  Paris  (1981),  123  pp. 

3.  A.  Bricaud,  A.  Morel,  and  L.  Prieur,  “Optical  Efficiency  Factors 
of  Some  Phytoplankters,"  Limnol.  Oceanogr.  28, 816  (1983). 

.4.  In  general  phytoplankton  scattering  will  not  satisfy  the  law 
Bc(\)  ~  X'n  for  a  constant  value  of  n  over  the  entire  visible 
spectrum.  This  law  is  used  here  only  to  provide  an  analytical 
representation  of  the  possible  spectral  behavior  over  this  very 
limited  portion  of  the  spectrum. 

15.  T.  J.  Petzold,  Volume  Scattering  Functions  for  Selected  Natu¬ 
ral  Waters,  Scripps  Institute  of  Oceanography,  Visibility  Lab¬ 
oratory,  San  Diego,  CA  92152,  SIO  Ref.  72-78  (1972). 


True  z0  (Maters) 

(b) 

Fig.  16.  Comparison  between  the  source  physical  depth  z0  re¬ 
trieved  by  inverting  Eq.  (36)  for  various  pigment  concentrations: 
(a)  /  *  50%;  (b)  f  =  90%. 

16.  H.  R.  Gordon,  “Ship  Perturbation  of  Irradiance  Measurements 
at  Sea.  1:  Monte  Carlo  Simulations,”  Appl.  Opt.  24,  4172 
(1985). 

17.  By  the  term  exact,  we  mean  in  the  sense  that  the  solutions  are 
given  in  terms  of  functions  which  can  be  evaluated  numerically, 
e  g.,  Chandrasekhar’s18  H  function. 

18.  S.  Chandrasekhar,  Radiative  Transfer  (Dover,  New  York, 
1960). 

19.  J.  P.  Elliott,  "Milne’s  Problems  with  a  Point  Source,”  Proc.  R. 
Soc.  London  Ser.  A  228, 424  (1955). 

20.  C.  Mark,  “Neutron  Density  Near  a  Plane  Surface,”  Phys.  Rev. 
72,658(1947). 

21.  R.  W.  Preisendorfer,  Hydrological  Optics.  Vol.  Ill:  Solu¬ 
tions.,  PB-259795/3ST  (National  Technical  Information  Ser¬ 
vice,  U.S.  Department  of  Commerce,  5285  Port  Royal  Rd., 
Springfield,  VA  22161  (1976). 

22.  R.  W.  Preisendorfer,  "On  the  Existence  of  Characteristic  Diffuse 
Light  in  Natural  Waters,"  J.  Mar.  Res.  18, 1  (1959). 

23.  M.  Abramowitz  and  I.  A.  Stegun,  Handbook  of  Mathematical 
Functions  (Dover,  New  York,  1965). 

24.  G.  Beardsley  and  J.  R.  V.  Zaneveld,  "Theoretical  Dependence  of 
the  Near- Asymptotic  Apparent  Optical  Properties  on  the  Inher¬ 
ent  Optical  Properties  of  Sea  Water,”  J.  Opt.  Soc.  Am.  59,  373 
(1969). 

25.  L.  Prieur,  “Transfertradiatif  dans  leseauxde  mer.  Application 
a  la  determination  de  parametres  optiques  caracterisant  leur 
teneur  en  substances  dissoutes  et  leur  contenu  en  particules,” 
D.Sci.  Thesis,  U.  Pierre  et  Marie  Curie  (1976),  243  pp. 


1148  APPLIED  OPTICS  /  Vol.  26.  No.  19  /  1  October  1987 


Appendix  2 


H.R.  Gordon,  Influence  of  Vertical  Stratification  on  the  Distribution  of  Irradiance  at  the 
Sea  Surface  from  a  Point  Source  in  the  Ocean,  Applied  Optics,  27,  2643-2645  (1988). 


Influence  of  vertical  stratification  on  the  distribution 
of  irradiance  at  the  sea  surface  from  a 
point  source  In  the  ocean 

Howard  R.  Gordon 

University  of  Miami,  Physics  Department,  Coral  Gables, 

Florida  33124. 

Received  25  November  1987. 

0003 -6935/88/ 1 32643 -03$02. 00/0. 

'©  1988  Optical  Society  of  America. 

In  this  Letter,  a  bio-optical1  model  describing  the  distribu¬ 
tion  of  irradiance  at  the  sea  surface  resulting  from  a  point 
source  embedded  in  a  homogeneous  ocean  is  extended  to  a 
stratified  ocean.  In  particular,  the  phytoplankton  pigment 
concentration-  C  is  allowed  to  vary  with  depth.  Briefly, 
following  Ref.  1,  the  ahsorption  a(  (z,A)  and  scattering  b((z,\) 
coefficients  of  the  phytoplankton  and  their  immediate  detri- 
tal  material  at  a  wavelength  A  are  given  by 

b,(z,X)  =  fi,  Wf(j)"K. 

a< (z.A)  =  0.06.4(  (A)C(z)1"'"'-’,  (1) 

where  z  is  depth,  A(  (A)  is  the  absorption  coefficient  of  plank¬ 
ton  pigments  normalized  to  that  at  440  nm,  and  Bc( A)  = 
0.3(550/A).  6(  and  afhave  units  of  m_l  when  A  is  in  nanome¬ 
ters  and  C  is  in  mg/m1.  The  beam  attenuation  coefficient 
c(z, A)  is  then 

c(z.A)  =  a„(A)  +  bJX)  +  a(  (X,z)  +  bc(\,z),  (2) 

where  au  (A)  and  bu( A)  are,  respectively,  the  scattering  coeffi¬ 
cients  for  pure  seawater.  The  scattering  phase  function  for 
the  particles  is  taken  from  Ref.  1  and  is  independent  of  depth 
and  pigment  concentration.  However,  the  total  scattering 
phase  function  (water  plus  particles)  does  depend  on  2  be¬ 
cause  the  relative  importance  of  scattering  by  particles  and 
water  depends  on  C(z).  The  radiative  transfer  equation  is 
solved  by  Monte  Carlo  techniques  to  describe  the  distribu¬ 
tion  of  irradiance  at  the  surface  produced  by  a  point  source  of 
unit  power  (1  VV/nm)  at  a  depth  z0.  Specifically,  the  irradi¬ 
ance  distribution  just  beneath  the  surface  Eu(R,z0)  is  deter¬ 
mined  as  a  function  of  R,  the  distance  measured  along  the 
surface  from  a  point  vertically  above  the  source.  Figure  1 
shows  the  results  of  four  Monte  Carlo  simulations  at  480  nm 
[A(  (480)  =  0.798]  in  which  the  source  was  at  20  =  30  m,  and 
the  pigment  concentration  was  (a)  0.5  mg/m3  from  0  to  10  m 
and  0.1  mg/m1  elsewhere,  (b)  0.5  mg/m3  from  10  to  20  m  and 
0.1  mg/m3  elsewhere,  (c)  0.5  mg/m3  from  20  to  30  m  and  0.1 
mg/m3  elsewhere,  and  (d)  a  uniform  distribution  of  pigments 
wi*h  a  concentration  equal  to  the  mean  concentration  in  (a), 

( b ) ,  and  (c)  above  the  source  (0.2333  mg/m3).  The  computa¬ 
tions  presented  in  Fig.  1  clearly  show  that  the  vertical  distri¬ 
bution  of  pigments  is  unimportant  in  determining  the  irradi¬ 
ance  distribution  at  the  surface;  only  the  mean  pigment 
concentration  and  the  source  depth  appear  to  be  relevant. 

In  Ref.  1  a  simple  analytical  model  was  presented  for 
computing  the  irradiance  distribution  from  a  point  source  of 
unit  power  embedded  in  a  homogeneous  ocean.  It  consisted 
of  computing  the  distribution  in  the  absence  of  scattering 
and  replacing  the  absorption  coefficient  by  an  effective  at¬ 
tenuation  coefficient  k,  i.e., 

EjR.z,,)  -  4^  ^exp(-AA),  (3) 

where  X2  =  R2  +  z[  and 

-tM/C.,  (4) 


with  /fsur f  and  K»  being  the  downwelling  irradiance  attenua¬ 
tion  coefficients,  respectively,  just  beneath  the  surface  and  in 
the  asymptotic  regime.  yUo)  is  given  by 

y(z0)  =  exp(-0.08  KI(ir,2„),  (5) 

K»urf  “  +  or  +  0.5fc*  +  (fc6),  bc,  (6) 

where  (bb)c  is  the  backscattering  probability  associated  with 
the  particle  scattering  phase  function.  K«  was  approximat¬ 
ed  by 

K-  *  K.uri  +  0.09c.  (7) 

The  model  reproduced  both  the  distribution  £u(/?,20)  in  the 
nondiffusion  regime  and  the  total  power  reaching  the  sur¬ 
face, 

f\ou.  “  j  2t Eu(R,z0)RdR  =  *  EJkzu), 

where  £2  is  the  exponential  integral,  with  high  accuracy. 
Since  Fig.  1  shows  that  Eu(R,z0)  is  apparently  independent  of 
the  stratification  in  C,  one  expects  that  the  analytical  model 
should  apply  equally  well  in  the  presence  of  the  stratification 
if  a  suitable  value  of  k  can  be  found.  Indeed,  the  following 
provides  a  straightforward  procedure  for  the  determination 
of  k  by  considering  an  equivalent  homogeneous  ocean. 

First,  an  effective  pigment  concentration  ( C)  is  defined  to 
be  the  (constant)  value  of  C,  which,  for  a  homogeneous  ocean 
and  given  source  depth  z0,  would  yield  the  same  value  for  the 
source  optical  depth  t  as  the  stratified  ocean,  i.e., 

r(A)  *»  £  c(z, X)dz  s  (c(A))z0,  (8) 

where  c(2,\)  is  found  by  inserting  the  pigment  profile  C(z) 
into  Eqs.  (1)  and  (2),  and  <c(A) >  is  the  beam  attenuation 
coefficient  for  the  equivalent  homogeneous  ocean.  Next 
c(2,X)  and  C(z)  in  Eqs.  (1)  and  (2)  are  replaced  by  <c(A)>  and 
(C> ,  respectively,  and  solved  for  (C>.  Finally,  this  value  of 
<C>  is  used  in  Eqs.  (3)-(7)  [which  also  employ  Eqs.  (1)  and 
(2)  with  C(z)  replaced  by  <  C)  ]  to  estimate  k?  An  example  of 
the  performance  of  the  analytical  model  using  this  procedure 
with  a  realistic  pigment  profile4  is  presented  in  Figs.  2  and  3. 
The  source  is  assumed  to  be  located  just  below  the  pigment 
maximum  at  z0  =  64  m  and  is  marked  with  a  dot  in  Fig.  2. 
The  results  of  a  Monte  Carlo  simulation  (noisy  curve)  and 
the  analytical  model  (smooth  curve)  for  Eu(R,z0 )  are  shown 
in  Fig.  3,  and  the  excellent  agreement  suggests  that  the 
analytical  model  can  provide  a  good  approximation  to  the 
time-independent  Green's  function  for  an  internal  source 
even  in  the  case  of  a  stratified  ocean. 

In  sum,  the  vertical  distribution  of  pigments  is  unimpor¬ 
tant  in  determining  the  horizontal  distribution  of  irradiance 
(and  P total)  from  a  point  source  in  the  water.  Also,  the  simple 
model  for  computing  the  irradiance  distribution  presented 
in  Ref.  1  provides  excellent  agreement  with  the  exact  results 
computed  using  Monte  Carlo  techniques  when  the  pigment 
concentration  is  replaced  by  (a  constant)  <C). 
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2.  By  the  term  pigment  concentration  we  mean  the  concentration 
(mg/m'1}  of  chlorophyll  a  and  all  chlorophyll  -like  pigmenta,  which 
absorb  in  the  same  spectral  bands  as  chlorophyll  a,  such  as 
phaeophytin  a,  and  which  are  contained  in  phytoplankton  or  in 
their  detrital  materials.  The  sum  of  the  concentrations  of  chloro¬ 
phyll  a  and  phaeophytin  a  is  frequently  used  as  an  indicator  of 
phytoplankton  biomass. 

3.  Note  that 

<0  *C*z~a 1  (*"  C{z)dz, 

Jo 

the  mean  pigment  concentration,  because  the  dependence  of 
b(  (z,\)  and  ai  (r,X)  on  C  in  Eq.  (1)  is  nonlinear. 

4.  E.  Swift,  U.  Rhode  Island;  personal  communication. 


Fig.  1.  Irradiance  distribution  at  the  surface  for  four  different 
pigment  profiles  having  the  same  total  amount  of  pigment  above  the 
source  (see  text  for  details). 


Pigment  Concentration  (mg/m^) 


Fig.  3.  Irradiance  distribution  at  the  surface  for  the  pigment  profile 
in  Fig.  2.  The  noisy  curve  represents  the  results  of  a  Monte  Carlo 
simulation,  and  the  smooth  curve  is  the  result  of  the  analytical  model 
(Eqs.  (3)-<7)J. 


Fig  2.  Pigment  concentration  as  a  function  of  depth  for  the  irradi¬ 
ance  distribution  in  Fig.  3.  The  dot  at  64  m  represents  the  position 
of  the  source. 
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a 

total  absorption  coefficient 

i 

m~l 

ap 

particle  absorption  coefficient 

m-1 

water  absorption  coefficient 

m-1 

b 

total  scattering  coefficient 

m-1 

bb 

total  backscattering  coefficient 

m"1 

h 

total  backscattering  probability 

particle  scattering  coefficient 

m~l 

bw 

water  scattering  coefficient 

m-1 

m 

total  volume  scattering  function 

m~1Ster~l 

/W) 

particle  volume  scattering  function 

m-1Ster-1 

/M©) 

water  volume  scattering  function 

m-1Ster~l 

c 

total  attenuation  coefficient  (a  -f  b) 

m-1 

cp 

particle  attenuation  coefficient  (ap  +  bp) 

m_1 

cw 

water  attenuation  coefficient  (a£  +  bw) 

m-1 

Ci 

attenuation  coefficient  of  component  i 

m_I 

< 

specific  attenuation  coefficient  c</C< 

m*mg-1 

C 

pigment  concentration 

mg  m~3 

Ci 

concentration  of  ith  constituent 

mg  m-3 

Do 

downwelling  distribution  function  (a>0  =  0) 

f 

direct  sun  fraction  of  £^(0) 

- 

F 

total  forward  scattering  probability  (1  -  6^) 

particle  forward  scattering  probability 

Fw 

water  forward  scattering  probability 

F0 

extraterrestrial  solar  irradiance 

mW  cm'Vm'1 

Ed(z) 

downwelling  irradiance  at  z 

mW  cm~2fxm~1 

Kd{z) 

attenuation  coefficient  for  Ed(z) 

m-1 

K 

Kd(  0) 

m"1 
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Kw 

Kp 

Kc 

Kx 
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(K)w 

(K)P 
( K)b 
(K)t 
A 
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P„(Q) 

<7J 

t 

T 

no 

TA 

TOl 

TR 

0 

t?0 

O/Q 

up 

ww 

z 

no 

Z*,  Zy 


pure  water  component  of  Kd( 0)  and  ^(z) 
water  component  of  K 
pigment  component  of  Kd(z) 
nonpigment-nonwater  component  of  Kd(z) 
mean  Kd(z)  from  surface  to  zx 0 
water  component  of  (K) 
particle  component  of  { K ) 

Lambert-Beer  value  of  ( K )  (( K)w  +  (K)p  ) 
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wavelength 

scattering  phase  function  (0/b) 
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Fresnel  transmittance  of  sea  surface 
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ozone  optical  thickness  of  atmosphere 

Rayleigh  optical  thickness  of  atmosphere 

scattering  angle 

solar  zenith  angle 

solar  zenith  angle  below  surface 

total  scattering  albedo  (6/c) 

particle  scattering  albedo  ( bp/cp ) 
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depth 

depth  at  10%  surface  irradiance 
surface  slope  components 


nm 
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m 
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Abstract 


The  transport  of  radiation  in  a  realistic  ocean-atmosphere  system  is  simulated,  and  the  results 
are  treated  as  experimental  data,  to  show  that  the  downwelling  irradiance  attenuation  coefficient 
just  beneath  the  surface  K  and  the  mean  irradiance  attenuation  coefficient  from  the  surface  to 
the  depth  where  the  irradiance  falls  to  10%  of  its  value  at  the  surface  ( K )  can  be  corrected  for 
the  geometric  structure  of  the  in-water  light  field  to  yield  quantities  that  are,  to  a  high  degree  of 
accuracy,  inherent  optical  properties.  These  geometry- corrected  quantities  are  shown  to  satisfy  the 
Lambert-Beer  law  to  a  reasonable  degree  of  accuracy,  with  the  largest  error  («  15%)  in  the  case  of 
( K )  arising  from  mixing  nonabsorbing  particles,  e.g.,  white  sand,  with  strongly  absorbing  water 
(wavelengths  >  600  nm).  This  near-validity  of  the  Lambert-Beer  law,  when  there  are  compelling 
reasons  to  believe  that  it  should  fail,  is  shown  to  result  from  three  independent  facts:  (1)  the 
dependence  of  the  diffuse  attenuation  coefficients  on  the  geometric  structure  of  the  light  field  can 
be  removed;  (2)  pure  sea  water  is  a  much  better  absorber  than  scatterer  at  optical  frequencies;  and 
(3)  the  phase  functions  for  particles  suspended  in  the  ocean  differs  significantly  from  that  of  pure 
sea  water.  Finally,  it  is  shown  that  extrapolation  of  the  corrected  diffuse  attenuation  coefficients  to 
the  limit  c  — *  cw  yields  quantities  which  are  within  2%  the  corresponding  quantities  that  would  be 
measured  for  an  ocean  consisting  of  pure  sea  water  with  the  sun  at  the  zenith  and  the  atmosphere 
removed. 


Introduction 


In  a  series  of  papers,  Smith  and  Baker  (Baker  and  Smith  1982,  and  references  therein)  have 
developed  a  “bio-optical”  model  for  relating  the  optical  properties  of  near-surface  ocean  water 
to  the  content  of  biological  material.  Specifically,  the  attenuation  coefficient  Kd  of  downwelling 
irradiance  Ed  defined  by  Kd  =  -(1/ Ej)dEd/dz,  where  z  is  depth,  is  related  to  the  phytoplankton 
pigment  concentration  C  through 


Kd  =  Kw  +  Kc{C)  +  Ka.  (1) 

C  is  the  concentration  (mg/m3)  of  chlorophyll  a  and  all  chlorophyll-like  pigments,  which  absorb 
in  the  same  spectral  bands  as  chlorophyll  a,  such  as  phaeophytin  a ,  and  which  are  contained  in 
phytoplankton  or  in  their  detrital  materials.  In  Equation  1,  the  Lambert-Beer  law  applied  to  Kd, 
Kw  is  the  contribution  to  Kd  from  the  water  itself,  Kx  is  the  contribution  from  material  suspended 
or  dissolved  in  the  water  and  not  covarying  with  C,  and  Kc(C)  represents  the  contribution  to 
Kd  from  phytoplankton  and  their  immediate  detrital  material.  This  decomposition  of  Kd  is  very 
attractive  for  the  optical  analysis  of  ocean  water  because  of  the  relative  ease  in  measuring  Ed, 
the  absence  of  the  requirement  for  absolute  radiometry  to  determine  Kd,  and  the  possibility  of 
measuring  Kd  remotely  (Austin  and  Petzold  1981,  Gordon  1982)  and  even  at  night  (Gordon  1987). 
To  utilize  it,  measurements  of  Kd  for  a  given  wavelength  and  from  a  variety  of  oceanic  waters  are 
plotted  as  a  function  of  C  and  the  minimum  envelope  of  the  resulting  curve  [(if ,*)„,(„]  is  assumed  to 
correspond  to  Kx  =  0.  Taking  the  limit  of  {Kd)mta  as  C  -►  0  yields  Kw.  Then,  Kc(C)  is  given  by 
( A'd)min  -  Kw.  (For  examples  of  this  procedure  see  Figure  1  in  Baker  and  Smith  (1982)  or  Figure  7 
in  Gordon  and  Morel  (1983).]  Assuming  that  Kc(C)  is  a  valid  for  all  waters,  Equation  1  can  then 
be  applied  to  specific  cases  to  estimate  Kx  from  Kd  and  C,  or  to  estimate  C  from  Kd  in  waters  for 
which  Kx  is  known  to  be  negligible.  These  latter  waters  are  usually  referred  to  as  “Case  1  waters,” 
and  are  defined  to  be  waters  for  which  the  optical  properties  are  controlled  by  phytoplankton  and 
their  immediate  detrital  material  (Gordon  and  Morel  1983,  Morel  and  Prieur  1977). 

Equation  1  has  been  criticized  by  Morel  and  Bricaud  (1981)  and  Stavn  (1988)  on  the  basis 
that,  unlike  the  absorption  coefficient  and  the  volume  scattering  function,  Kd  is  not  an  inherent 
optical  property  of  the  medium  Preisendorfer  (1961).  This  is  because  it  depends  on  the  depth  and 
on  the  geometric  structure  of  the  in-water  light  field,  as  well  as  on  the  inherent  optical  properties 
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of  the  medium.  Since  a  given  Kj  is  unique  only  to  the  particular  situation  in  which  it  is  measured, 
and  there  is  no  reason  to  expect  that  the  three  components  of  Kj  will  vary  in  the  same  manner 
with  depth  and  with  the  structure  of  the  light  field,  it  is  correctly  asserted  that  Eq.  1  is  only 
an  approximation.  However,  Gordon,  Brown  and  Jacobs  (1975)  have  shown  with  Monte  Carlo 
simulations  of  the  in-water  light  field,  that  for  simple  modes  of  illumination,  i.e.,  a  sky  of  uniform 
radiance  or  a  parallel  beam  of  irradiance  incident  at  an  angle  i?o  with  the  vertical,  the  dependence 
of  Kd  on  the  structure  of  the  light  field  can  be  removed  v/ithout  any  knowledge  of  the  optical 
properties  of  the  medium.  Because  of  this,  they  called  the  light  field-corrected  Kd  a  quasi-inherent 
optical  property  of  the  medium.  Later,  Gordon  (1976)  showed  that  the  correction  factor  required 
to  remove  the  light  field  dependence  from  Kd  could  be  reasonably  accurately  computed  by  knowing 
only  the  relative  amounts  of  skylight  and  direct  sunlight  incident  on  the  sea  surface  in  the  spectral 
band  in  question.  Later,  Baker  and  Smith  (1979)  directly  verified  that  for  turbid  water  under 
clear  skies  with  the  sun  near  the  zenith,  Kd  was  nearly  independent  of  i?0.  Finally,  Gordon  (1980) 
demonstrated  that  for  a  stratified  ocean  with  sun  at  the  zenith,  the  value  of  Kd  at  a  given  depth  z 
depended  mostly  on  the  inherent  optical  properties  of  the  medium  at  that  depth,  i.e.,  Kd  is  a  local 
property  of  the  medium.  These  observations  concerning  the  quasi-inherent  nature  of  Kd  suggest 
that  if  Kd- values  corrected  for  variations  in  the  illumination,  e.g.,  corrected  so  that  the  resulting 
Kd  is  that  value  which  would  be  measured  in  the  same  ocean  illuminated  by  the  sun  at  the  zenith 
in  the  absence  of  the  atmosphere  (no  sky  light),  are  used  in  Eq.  1  the  error  resulting  from  the  fact 
that  Kd  is  not  a  true  inherent  optical  property  would  be  considerably  reduced.  However,  a  residual 
error  would  still  remain  in  Eq.  1  because  Kd  depends  on  depth. 

In  this  paper  the  earlier  computations  (Gordon,  Brown  and  Jacobs  1975)  of  Kd  are  extended 
to  cases  of  more  realistic  illumination  of  the  surface  and  a  more  realistic  model  of  the  inherent 
optical  properties  of  the  water  itself.  The  results  confirm  that  when  Kd  for  the  water  is  measured 
either  just  beneath  the  surface,  or  an  average  value  is  determined  between  the  surface  and  the 
depth  where  the  surface  irradiance  is  reduced  to  10%  of  its  value  at  the  surface,  the  resulting  Kd’ s 
can  be  corrected  to  yield  a  quantity  that  can  be  directly  expressed  in  terms  of  the  inherent  optical 
properties  of  the  medium.  Furthermore,  this  quantity  is  shown  to  depend  nearly  linearly  on  the 
inherent  optical  properties  of  the  medium,  which  are  linearly  additive  over  the  constituents,  and  so 
it  satisfies  the  Lambert-Beer  law  with  reasonable  accuracy.  Finally,  the  results  show  that  the  value 
of  Am  determined  by  extrapolation  of  Kd  to  C  =  0  in  Case  1  waters  is  in  fact  very  nearly  equal  to 
the  value  of  Kd  that  would  be  measured  in  an  ocean  consisting  of  only  pure  sea  water. 
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Optical  Model  of  the  Ocean- Atmosphere  System 

Our  study  of  the  efficacy  of  Equation  1  is  based  on  simulating  the  transfer  of  radiation  in 
the  ocean-atmosphere  system  by  Monte  Carlo  techniques.  Such  simulations  provide  the  distribu¬ 
tion  of  radiation  in  the  entire  system.  The  radiation  field  is  treated  as  experimental  data,  albeit 
data  collected  under  carefully  controlled  conditions,  a  cloud  free  sky  and  a  homogeneous  ocean  of 
precisely  known  inherent  optical  properties,  from  which  the  irradiances,  Kd,  etc.,  can  be  derived 
as  a  function  of  the  inherent  optical  properties  of  the  ocean  and  ultimately  as  a  function  of  the 
constituent  concentrations.  To  accomplish  this,  optical  models  of  the  ocean  and  the  atmosphere 
that  are  as  realistic  as  possible  are  required.  Such  models  are  described  below.. 

A.  The  Ocean 

Water  find  its  constituents  influence  Kd  through  their  effect  on  the  inherent  optical  properties  of 
the  medium:  a  the  absorption  coefficient;  and  /3(0)  the  volume  scattering  function.  For  simplicity, 
we  limit  the  modeling  of  these  properties  to  Case  1  waters.  Moreover,  we  also  limit  the  model 
to  those  waters  for  which  optically  active  dissolved  oxganic  materials  —  yellow  substances  —  are 
absent.  The  rational  for  this  is  that  the  effect  of  yellow  substance  absorption  on  the  results  would 
be  identical  to  that  of  increasing  the  absorption  coefficient  of  pure  sea  water  by  an  appropriate 
amount.  Therefore,  the  medium  is  described  by  an  absorption  coefficient  given  by 

a  =  aw  +  ap,  (2) 

where  the  subscripts  w  and  p  here,  and  hereafter,  refer  to  the  contribution  from  water  and  suspended 
particles,  respectively,  and  a  volume  scattering  function  given  by 

0(0)  =  &.(©) +  &(©). 

where  0  is  the  scattering  angle.  The  total  scattering  coefficient  b  is  defined  by 

b  =  2ir  f  /3(0)sin0d0, 

Jo 

and  b  =  bw  +  bp.  The  total  attenuation  coefficient  is  defined  to  be 

c  =  Cu,  -f-  cp  =  a  -f  6.  (5) 


(3) 

(4) 
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It  is  convenient  to  define  two  other  auxiliary  parameters,  the  single  scattering  albedo  =  b/c 
and  the  scattering  phase  function  P(0)  =  @{0)fb.  The  value  of  u>o  is  the  probability  that  at  each 
interaction  within  the  medium  the  photon  will  be  scattered  rather  than  absorbed,  while 

,&3 

2jt  /  P(0)sin0d0 

JQi 

is  the  probability  that,  when  a  photon  is  scattered,  the  scattering  angle  will  be  between  0j  and 
02-  A  quantity  which  is  needed  later,  the  backscattering  probability  bb,  is  defined  by 

bb  =  2tt  f  P(0)sin0d0,  (6) 

J  w/2 

and  is,  therefore,  seen  to  be  the  probability  that  a  photon  is  scattered  through  an  angle  larger  than 
90°  .  The  backscattering  coefficient,  bb,  is  then  given  by  bb  =  bbb. 


It  is  easy  to  verify  from  the  radiative  transfer  equation  that  for  a  homogeneous  ocean  the 
internal  radiation  field  is  completely  specified  by  providing  c,  u>o,  P(0),  and  the  distribution  of 
radiation  incident  on  the  sea  surface,  and  also  that  the  field  depends  on  depth  only  through  the 
product  cz.  P(0)  and  oj0  can  be  related  to  the  similar  water  and  particle  quantities  PW(Q),  Pp(0), 
u>w,  and  u>p,  through 

~  cp/cT+  i"’’  (7) 


UJqP  _  up^>p(cp/c'°)  d*  uv)Py> 

cp/  41 1 

Thus,  in  this  simple  model,  given  cw,  cp,  <mw,  wp,  Pw(©)»  and  Pp(0),  representing  the  inherent 
optical  properties  of  the  water  and  the  pair  tides,  u>o  and  P(0)  are  specified  by  the  ratio  cp/cw. 
This  ratio  is  proportional  to  the  particle  concentration. 


Experimental  measurements  must  be  used  to  provide  a  realistic  parameterization  of  the  optical 
properties  ww,  ojp,  Pw(0),  and  Pp(0).  The  model  used  here  is  identical  to  that  developed  by  the 
author  (Gordon  1987)  to  study  the  propagation  of  irradiance  from  a  point  source  embedded  in  the 
ocean.  Briefly,  the  absorption  coefficient  aw  has  been  inferred  from  measurements  of  downwelling 
and  upwelling  irradiance  in  oligotrophic  waters  such  as  the  Sargasso  Sea  (Baker  and  Smith  1982, 
Morel  and  Prieur  1977,  Prieur  and  Sathyendranath  1981),  and  the  scattering  coefficient  bw  and 
the  volume  scattering  function  0W(0)  have  been  measured  directly  for  pure  water  and  for  saline 
solutions  of  pure  water  corresponding  to  salinities  between  35  and  39  ppt  by  Morel  (1974).  The 
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resulting  aw,  bw  and  u>w  are  given  in  Table  1  for  the  wavelengths  used  in  the  present  computations. 
Note  that  these  values  of  u>w  represent  the  upper  limit  of  the  scattering  albedo  for  water  plus  any 
dissolved  material  such  as  yellow  substances,  since  the  small  concentrations  of  dissolved  material 
typically  found  in  sea  water  can  contribute  to  the  absorption  coefficient  but  not  to  the  scattering 
coefficient.  The  scattering  phase  function  for  pure  sea  water  is  taken  from  Morel  (1974). 


The  optical  properties  of  the  suspended  particles  for  Case  1  waters,  can  be  related  to  the 
pigment  concentration.  The  scattering  coefficient  of  particles  at  550  nm,  5P(550),  is  nonlinearly 
related  to  the  pigment  concentration  C  through  (Morel  1980) 

bp  =  BcCon,  (9) 


where  6P(550)  is  in  m_1  and  C  is  in  mg/m3,  (See,  also,  Gordon  and  Morel  (1983)).  The  constant 
Be,  the  scattering  coefficient  at  a  pigment  concentration  of  1  mg/m3,  ranges  from  0.12  to  0.45  and 
has  an  average  value  of  0.30.  The  variation  in  Be  is  due  to  the  natural  variability  of  scattering  over 
the  various  species  of  phytoplankton,  as  well  as  a  variability  in  scattering  by  the  detrital  particles 
associated  with  the  phytoplankton.  Similarly  the  absorption  coefficient  of  the  particles  has  been 
studied  as  a  function  of  C  by  Prieur  and  Sathyendranath  (1981)  yielding  for  C  <  10  mg/m3: 

op(A)  =  0.06^c(A)C'0-602,  (10) 


where  ap(A)  is  in  m  1  and  C  is  in  mg/m3.  In  this  equation  .Ac (A)  is  the  absorption  coefficient  of 
phytoplankton  normalized  to  440  nm,  t.e., 


Ac(  A)  = 


ap(A) 

ap(440)’ 


(11) 


These  nonlinear  relationships  between  bp  and  C  and  ap  and  C  are  believed  to  be  due  to  a  systematic 
variation  in  the  ratio  of  the  concentration  of  phytoplankton  to  that  of  detrital  material  as  a  function 
of  the  concentration  of  phytoplankton  (Hobson,  Menzel  and  Barber  1973,  Smith  and  Baker  1978a). 
The  relative  absorption  of  phytoplankton  ^c(A)  deduced  by  Prieur  and  Sathyendranath  (1981) 
agrees  well  with  absorption  measurements  made  on  phytoplankton  cultures  by  Sathyendranath 
(1981).  Note  that  U|>(A)  includes  both  phytoplankton  and  their  detrital  material  and  thus  represents 
the  absorption  of  all  components  other  than  the  water  itself. 
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Since  bp( A)  and  ap(A)  vary  with  pigment  concentration  in  nearly  the  same  manner,  6p(A)/ap(A) 
is  nearly  independent  of  the  pigment  concentration,  i.e., 


ap(A)  ~ 


16.6 


Bc{  A) 

Ac(\y 


(12) 


This  provides  an  estimate  of  wp(A),  and  shows  that  this  quantity  is,  in  the  first  approximation, 
independent  of  the  pigment  concentration.  At  550  nm,  where  Be  is  known,  this  yields  u>p(55o)  = 
0.933  in  good  agreement  with  the  range  for  those  measured  by  Bricaud,  Morel  and  Preiur  (1983) 
for  four  species  of  cultured  phytoplankton:  0.89  <  u>p(550)  <  0.97.  In  order  to  fix  reasonable  values 
of  u)p( A)  at  the  other  wavelengths  of  interest,  the  variation  of  Be  with  A  is  required.  Following 
Gordon  (1987)  we  assume  Bc(A)  obeys  a  power  law  with  wavelength,  i.e.,  Bc( A)  oc  A-n,  and 
take  n  =  +1.  This  yields  Bc( 480)  ss  0.34  and  Be (440)  ss  0.38.  The  resulting  values  of  wp(A)  and 
w«u(A)  used  in  the  computations  are  provided  in  Table  2.  It  should  be  noted  that  the  assumption 
Bc( A)  oc  A-1  often  overestimates  the  dependence  of  be  on  A  since  the  scattering  by  absorbing 
particles,  e.g.,  phytoplankton,  tends  to  be  depressed  in  the  pigment  absorption  bands  (Bricaud, 
Morel  and  Preiur  1983).  This  depression  of  scattering  would  make  h/p(440)  and  u/p(480)  in  smaller 
than  given  in  Table  2;  however,  the  effect  is  not  large,  e.g.,  changing  n  from  +1  to  —1  only  reduces 
wp(440)  from  0.86  to  0.80.  To  insure  that  wide  departures  of  up  from  those  used  in  Table  2  do  not 
influence  the  results  of  this  work,  simulations  also  have  been  carried  out  for  o»p(480)  =  0.5,  0.7,  and 
0.99. 


The  particle  phase  function  is  the  most  difficult  quantity  to  parameterize  because  it  requires 
the  individual  phase  functions  of  the  plankton  and  the  detrital  material,  neither  of  which  have 
ever  been  measured  in  the  field.  Thus,  we  must  rely  on  measurements  of  the  total  particle  phase 
function  (plankton  plus  detrital  material).  Measurements  of  the  volume  scattering  function  at  530 
nm  have  been  made  for  waters  in  several  locations  with  very  different  turbidities  (total  scattering 
coefficients)  by  Petzold  (1972).  When  the  scattering  by  pure  sea  water  is  subtracted,  the  resulting 
particle  phase  functions  are  very  similar,  having  a  standard  deviation  which  is  within  about  30% 
of  the  mean,  over  waters  for  which  the  particle  scattering  coefficient  varied  over  a  factor  of  50. 
This  mean  particle  phase  function  derived  from  Petzold’s  measurements  is  adopted  for  this  study 
and  designated  by  the  symbol  “M.n  Also,  two  other  particle  phase  functions  are  used  to  represent 
the  extremes  of  the  phase  functions  given  by  Petzold’s  measurements.  These  are  (1)  the  mean  of 
three  phase  functions  measured  in  the  turbid  waters  of  San  Diego  Harbor  and  designated  by  “T,” 
and  (2)  a  phase  function  measured  in  the  clear  waters  of  the  Tongue  of  the  Ocean,  Bahamas,  and 
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designated  by  “C."  The  three  particle  phase  functions  Me  shown  in  Figure  1  along  with  the  phase 
function  for  scattering  by  the  water  itself  (Pw(0)).  It  is  seen  that  these  model  phase  functions 
differ  principally  in  their  scattering  at  angles  greater  than  25° .  The  backscattering  probabilities, 
(bt,)p,  associated  with  them  are  0.0120,  0.0144,  and  0.0181,  respectively  for  C,  M,  and  T. 

This  completes  the  specification  of  the  quantities  needed  for  the  simulation:  u„,  u>p,  PW(Q), 
ar.d  Pp(0).  Varying  the  parameter  cp/cw  from  0  to  oo  results  in  models  v'hicli  range  from  a  particle 
free  ocean  to  an  ocean  in  which  the  optical  properties  of  the  particles  are  completely  dominant.  This 
parameter  can  be  related  to  the  pigment  concentration  through  the  bio-optical  model  by  noting 
that  cp  =  ap  - f  bp  and  using  Equations  9  and  10.  The  result  is 

^T~7(A)C»«,  (13) 

where  7(A)  =  22.4,  18.5,  and  4.9  at  440,  480,  and  550  nm,  respectively. 


B.  The  Atmosphere 


The  atmosphere  influences  K4  by  distributing  a  portion  of  the  near-parallel  solar  beam  over 
the  entire  upward  hemisphere,  i.e.,  in  producing  sky  light  from  direct  sunlight.  In  order  to  simulate 
the  angular  distribution  of  radiation  entering  the  ocean,  an  atmospheric  model  is  required.  This 
model  atmosphere  consisted  of  fifty  layers  and  included  the  effects  of  aerosols,  ozone,  and  Rayleigh 
scattering,  vertically  distributed  according  to  data  taken  from  the  work  of  Elterman  (1968).  The 
aerosol  phase  functions  were  computed  by  Fraser  (R.  Fraser,  NASA/GSFC,  Personal  Communica¬ 
tion)  from  Mie  theory  using  the  Deirmendjian  (1969)  Haze  C  size  distribution.  This  model  simulates 
optical  properties  of  the  cloud- free  atmosphere  only. 


Computations  and  Properties  of  K& 

The  model  presented  above  is  used  to  specify  hypothetical  ocean- atmosphere  systems  for  which 
the  radiative  transfer  equation  is  solved  by  Monte  Carlo  techniques  to  provide  the  internal  radiation 
field.  The  radiation  field  is  then  treated  as  experimental  data  and  used  to  derive  the  quantities 
of  interest.  In  most  of  the  simulations  the  ocean  is  flat.  The  various  model  oceans  are  specified 
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by  allowing  the  pigment  concentration  C  to  vary  from  0  to  about  4.5  mg/ m3  which  in  turn  causes 
Cp/cw  to  vary  according  to  Eq.  13.  This  variation  in  cp/cw  then  induces  variations  in  u0  and  P(0) 
determined  by  Equations  7  and  8.  Figure  2  shows,  for  example,  the  change  in  the  shape  of  the  total 
phase  function  at  480  nm  as  cp/cw  is  varied  from  0  to  100.  Note  how  the  phase  function  deviates 
strongly  from  that  of  pure  water  (Rayleigh  scattering)  even  when  cp/cw  =  1,  i.e.,  even  when  the 
total  attenuation  is  shared  equally  between  water  and  particles. 

From  the  Monte  Carlo  solution  of  the  transfer  equation  the  irradiance  as  a  function  of  depth 
can  be  estimated.  Actually,  Ed(r),  where  r  is  called  the  optical  depth  (r  =  cz),  is  computed  in  the 
simulations.  The  irradiance  attenuation  coefficient  Kd  can  then  he  determined  from 

£i  =  'N*W>1  (14) 

c  dr 

by  numerical  differentiation.  Individual  values  of  Kd  computed  at  the  midpoint  of  the  euphotic 
zone  using  particle  phase  function  WT”  and  Cp/cu,  very  large,  but  in  the  absence  of  the  atmosphere, 
agree  well  with  the  values  computed  by  Kirk  (1984).  Figure  3  provides  some  samples  of  the 
resulting  profiles  of  the  irradiance  attenuation  coefficient.  In  these  examples  ifj/c  is  computed  for 
cP/cw  =  0,  1.4,  and  5.6  for  particle  phase  function  “M”  at  440  nm  with  the  sun  at  the  zenith.  The 
deepest  computed  point  for  each  profile  corresponds  to  r  =  9.  These  particular  profiles  represent 
a  reasonably  clear  ocean,  i.e.,  C  <  0.1  mg/ms.  The  immediate  conclusion  to  be  drawn  from  these 
simulations  is,  as  mentioned  in  the  introduction,  that  Kd  is  dependent  on  the  depth  even  for  a 
homogeneous  ocean.  Also,  Kd  increases  more  rapidly  with  z  at  larger  values  of  cp/cw.  In  fact, 
from  the  surface  to  z  —  100  m,  Kd  increases  by  2.5%,  10%  and  20%  for  Cp/cw  =s  0,  1.4,  and  5.6, 
respectively.  It  is  also  seen  that  the  rate  of  increase  in  Kd  is  more  rapid  near  the  surface.  Analysis 
of  Kd  just  beneath  the  surface  shows  that  the  rate  of  change  of  Kd  with  depth,  i.e.,  d^Kd*/e^ ,  is 
approximately  (0.02  ±  0.005)  X  c. 

These  computations  confirm  the  argument  that  Kd  cannot  be  considered  an  inherent  optical 
property  because  it  depends  on  depth,  and  also  the  suspicion  that  the  manner  of  the  dependence 
of  Kd  on  depth  is  a  function  of  the  particle  concentration.  [An  exception  to  this  of  course  is  the 
asymptotic  light  field  (z  — »  oo)  for  which  it  has  been  shown  (Preisendorfer  1959)  that  Kd  becomes 
independent  of  depth.]  Thus,  if  one  attempts  to  use  Kd  as  an  inherent  optical  property,  it  is 
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necessary  to  specify  in  some  manner  the  depth  at  which  the  measured  value  applies.  In  the  present 
work,  we  focus  on  the  irradiance  attenuation  coefficient  (K)  just  beneath  the  surface,  i.e., 

K  =  lim  Kd(r\  (15) 

T— »0 

and  on  the  average  diffuse  attenuation  coefficient  ((K))  over  the  upper  half  of  the  euphotic  zone, 

(K)  _  ln(Ed{Tio)/ Ed( 0))  ^16\ 

C  Tio 

where  t10  is  the  optical  depth  for  which  Ed  falls  to  10%  of  its  value  just  beneath  the  surface 
[Ed(Tl0)/Ed(  0)  =  0.1]. 

In  order  to  simulate  all  cloud  free  situations,  the  computations  have  been  carried  out  for  solar 
zenith  angles  of  0,  20°  ,  25° ,  30°  ,  40°  ,  60°  and  80° .  Also,  to  simulate  a  totally  overcast  sky, 
each  ocean  model  has  been  studied  with  the  atmosphere  removed  and  a  totally  diffuse  light  field 
incident  on  the  sea  surface.  Thus,  only  situations  with  broken  clouds  are  not  considered  in  this 
work.  A  sky  with  broken  clouds  is  particularly  difficult  to  examine  because  the  ladiation  field  is 
no  longer  independent  of  the  observer’s  horizontal  position  in  the  medium. 

Figures  4  and  5  provide  the  computations  of  K/c  and*{jRT)/c ,  respectively,  as  a  function  of  u>0  for 
334  simulations  comprising  a  variety  of  t?o ’s  and  C's  for  each  particle  phase  function  and  wavelength. 
Note  that  1  —  u/o  =  afc,  so  these  figures  relate  K  and  (K)  to  the  absorption  coefficient  a.  Based 
on  the  number  of  photons  contributing  to  K  and  ( K )  ,  the  statistical  uncertainty  6K/c  in  jRf/c  is 
a  ±0.017,  while  the  relative  error  6{K)  in  {K)  is  approximately  ±0.006,  i.e.,  8{K)  /( K )  »  ±0.006. 
Although  a  strong  trend  of  increasing  K/c  and  (K)/c  with  an  increasing  absorption  component  in 
the  total  attenuation  is  observed,  it  is  clear,  as  expected,  that  the  variation  in  K  and  {K)  cannot 
be  explained  solely  on  the  basis  of  the  total  absorption  and  scattering  coefficients  of  the  medium 
alone. 

In  order  to  proceed  further  it  is  useful  to  review  the  results  of  earlier  investigations.  Gordon, 
Brown  and  Jacobs  (1975)  found  that  the  dependence  of  Kd(r)  on  the  scattering  phase  function 
could  be  approximately  removed  by  expressing  Kd{r)  as  a  function  of  usqF,  where  F  is  the  forward 
scattering  probability  (F  =  1  —  6*),  rather  than  u;0  alone.  Also,  they  found  that  the  effect  of  the 
nature  of  the  illumination  of  the  ocean  on  Kd(r)  could  be  understood  by  examining  Kd{r)/ Dq(t), 
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where  Z?o(t)  was  the  downwelling  distribution  function  (Preisendorfer  1961)  for  a  totally  absorbing 
ocean  with  the  same  surface  illumination,  i.e., 


=  tw  <17) 

with  u>o  =  0,  where  Eod  is  the  downwelling  scalar  irradiance.  In  the  revised  suggested  notation 
(Morel  and  Smith  1982)  for  optical  oceanography  Do(t)  =  1/J ld  for  u>o  =  0,  where  JId  is  the 
“average  cosine”  of  the  downwelling  light  field  evaluated  just  beneath  the  surface.  In  the  study 
in  Gordon,  Brown  and  Jacobs  (1975)  there  was  no  atmosphere  over  the  ocean,  and  in  that  case 
Do(t)  =  1/  cos  d0uM  where  t?ou»  is  the  solar  zenith  angle  measured  beneath  the  sea  surface.  In 
the  present  simulations  Do  depends  on  wavelength  because  the  amount  of  skylight  produced  by 
scattering  in  the  atmosphere  is  a  function  of  wavelength.  Also,  the  relative  amounts  of  skylight 
and  direct  sunlight  depends  on  t?o-  Therefore,  Dq  has  been  computed  at  each  wavelength  and 
for  each  solar  zenith  angle  by  directly  solving  the  transfer  equation  for  the  given  A  and  $o  with 
u>0  =  0.  The  results  of  this  computation  for  Do  just  beneath  the  surface  (t  =  0)  are  provided  in 
Table  3.  Note  that,  as  expected,  D0  usually  increases  with  increasing  t?0;  however,  for  440  nm  the 
contribution  from  the  increasing  amount  of  skylight  compared  to  direct  sunlight  from  t?0  =  60°  to 
t?o  =  80°  actually  causes  a  small  decrease  in  Do-  Also,  for  t?o  <  60°  the  difference  between  I?o(0) 
and  1/  cosi?o«>  is  usually  less  than  about  3%.  Dq{t)  also  depends  on  r;  however,  this  dependence 
is  of  little  interest  here. 

Applying  the  observations  from  previous  studies  to  the  computations  in  Figures  4  and  5, 
Figures  6  and  7  provide  K/cDq  and  (K)/cDo  ,  respectively,  as  a  function  of  1  —  ujoF,  where  Do 
is  the  value  of  -Do(0)  taken  from  Table  3.  It  is  seen  that  when  the  computations  are  presented  in 
this  manner,  K/cDo  and  (K)/cDo  fall  on  what  appear  to  be  universal  curves.  The  curves  on  the 
figures  are  least-squares  fits  of  the  points  to 

K  2 

7d~o  =  S  *»  ~  U,‘>-P)n*  (18) 

and 

(19) 

n=l 

with  k\  —  1.0617,  &2  =  —0.0370,  ( k)\  —  1.3197  =  —0.7559,  and  (k) 3  =  0.4655.  The  average 

error  in  the  least-squares  fit  to  Equations  18  and  19,  is  1.8%  and  2.2%,  respectively.  (Replacing 
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jPo(O)  by  X>o(Tio)  in  Equation  19  provides  no  significant  increase  in  the  quality  of  the  expansion.) 
Also,  a  linear  fit  of  K/cD0  to  (1  -  waF)  is  almost  as  good  as  Eq.  18,  i.e., 

=  1.0395(1  -  u>oF),  (20) 

cDq 

which  implies  that 

-£-  =  1.0395 (a -Ms),  (21) 

JJq 

with  an  average  error  of  2.5%.  It  is  seen  that  the  points  on  Figures  6  and  7  with  1  -  wqF  >  0.85 
do  not  fit  Equations  18  and  19  quite  as  well  as  the  rest.  This  steins  from  the  fact  that  these 
points  correspond  to  pure  sea  water,  the  phase  function  of  which  differs  considerably  from  an  ocean 
containing  particles  (see  Figure  2).  The  K/cD0  -  (1  -  u>oF)  relationship  computed  for  an  ocean 
free  of  particles  is  presented  in  Figure  8  (for  t?o  =  0)  and  is  seen  to  differ  considerably  from  that  in 
Figures  6  and  7.  Since  the  minimum  value  of  (1  —  mu>Fw)  is  0.85  (near  400  nm),  and  over  the  range 
0.85  <  (1  -  u >WFW)  <  1  the  K/cD0  -  (1  -  u0F)  relationship  for  water  and  for  the  strongly  forward 
scattering  particles  are  very  similar,  the  computations  for  the  model  ocean  all  fall  very  near  the 
universal  curves  even  though  there  is  a  large  variation  in  the  shape  of  the  scattering  phase  function. 

The  above  analysis  shows  that  K/Dq  and  ( K)/Do  can  be  written  as  explicit  algebraic  functions 
of  the  inherent  optical  properties  c,  ujq  and  F  (independently  of  the  geometrical  nature  of  the 
light  field)  with  an  accuracy  that  is  likely  better  than  the  accuracy  with  which  K  or  ( K )  can  be 
measured.  Therefore  we  are  justified  in  regarding  the  quantities  K/Do  and  ( K)/D0  as  inherent 
optical  properties.  Physically,  K/Dq  and  ( K)/Dq  are  the  values  of  K  and  (K)  that  would  be 
measured  for  a  given  ocean  if  the  atmosphere  were  removed,  the  sun  placed  at  the  zenith,  and  the 
sea  surface  rendered  absolutely  flat.  In  such  a  well-defined  setting,  it  should  not  be  surprising  that 
K  and  (K)  can  be  considered  to  be  inherent  optical  properties.  What  is  surprising,  however,  is 
that  the  results  of  measurements  in  reed  situations  can  be  transformed  to  this  ideal  setting  through 
the  simple  division  by  Do . 

To  consider  applying  this  result  to  a  real  ocean,  the  effect  of  surface  roughness  on  this  simple 
observation  must  be  examined.  In  order  to  include  surface  waves  in  the  radiative  transfer  code  a 
statistical  model  of  the  waves  is  required.  For  simplicity,  we  assume  that  the  surface  roughness 
has  no  preferred  direction,  i.e.,  the  structure  of  the  surface  is  independent  of  the  wind  direction. 
Then  using  the  measurements  of  Cox  and  Munk  (1954)  the  probability  density  that  the  sea  surface 
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at  a  given  point  has  slope  components  zx  and  zv,  respectively,  in  the  x  and  y  directions  is  given 
approximately  by 

where  <r2 ,  the  slope  variance,  is  related  to  the  wind  speed  V  (in  m/s)  through 


a1  =  0.003  + 0.00512  V. 


The  rough  surface  described  by  p(zx,  zy)  is  incorporated  into  the  Monte  Carlo  radiative  transfer 
code  used  in  this  work  in  the  manner  similar  to  that  described  by  Plass,  Kattawar  and  Guinn 
(1975).  A  complete  examination  of  the  effect  of  surface  roughness  on  K  and  ( K )  requires  a  signif¬ 
icant  computational  effort;  however,  only  a  few  computations  are  required  to  show  that  the  basic 
result  above  —  division  by  Do  renders  K  and  (K)  inherent  optical  properties  —  is  still  valid  for 
an  ocean  with  waves.  A  sample  of  the  computations  carried  out  is  presented  in  Table  4  which 
provides  computations  of  D0,  K  and  (K)  as  a  function  of  the  surface  roughness  at  480  nm  for 
t?o  =  60°  and  cvjcw  —  12.3  (C  «  0.5  mg/m3).  Note  the  slow  increase  in  Do  with  a  indicating 
an  increasingly  diffuse  incident  light  field  beneath  the  surface  as  the  roughness  increases.  This 
increases  K  with  increasing  roughness;  however,  division  of  K  by  Do  provides  a  quantity  that  is 
nearly  independent  of  the  surface  roughness.  Interestingly,  the  effect  of  surface  roughness  on  both 
(K)/c  and  {K)/cD0  is  small  (<  3%)  up  to  wind  speeds  of  17  m/s.  These  computations  suggest 
that  K / Do  and  (K)/ D0  remain  inherent  optical  properties  even  in  the  presence  of  surface  waves; 
however,  the  value  of  Dq  used  to  form  these  ratios  must  be  that  which  is  valid  in  the  presence  of 
the  rough  surface. 


Determination  of  Do  from  field  measurements  requires  the  radiance  distribution  incident  on 
the  sea  surface.  This  can  be  quantitatively  determined  using  a  camera  equipped  with  a  fisheye  lens 
(Smith  1974,  Smith,  Austin  and  Tyler  1970);  however,  analysis  of  the  resulting  sky  photographs  is 
not  simple.  Earlier,  (Gordon  1976)  I  proposed  a  simple  scheme  for  estimating  Do.  Briefly,  if  Ed(i) 
is  the  irradiance  incident  on  the  sea  surface  from  source  i;  e.g.,  direct  sunlight,  skylight,  clouds, 
etc.,  then  it  is  easy  to  show  that 


Si  A >(i)t(i)Ed(i) 
Si  *(0^(0  ’ 


(22) 
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where  t(i)  is  the  irradiance  transmittance  for  light  from  source  *  and  Do(t )  is  the  value  of  Do  that 
would  result  from  source  i  acting  alone.  For  a  cloud  free  atmosphere  the  only  sources  are  the  sun 
and  the  sky  and  Equation  22  reduces  to 

Do  =  /Do(Sun)  +  (l-/)D0(Sky),  (23) 

where  /  is  the  fraction  of  direct  sunlight  in  the  incident  irradiance  transmitted  through  the  interface, 
i.e., 

_ f(Sun)  Ed(  Sun) _ 

t(Sun)  Ej(Sun)  +  i(Sky)Ed(Sky)' 

If  skylight  is  assumed  to  have  a  uniform  radiance  distribution,  i.e.,  radiance  (brightness)  indepen¬ 
dent  of  direction  of  viewing,  then  Equation  23  simplifies  to 

Do  =  ~~  +  U969(l  -  /).  (24) 

COS  1/0  uf 

Given  i?0  the  only  unknown  in  Equation  24  is  /.  This  can  be  estimated  by  placing  an  irradiance 
meter  above  the  surface,  measuring  the  total  incident  irradiance  Ed(Sun)  -f  Ed(Sky),  and  then 
measuring  the  sky  irradiance,  Ed(Sky),  by  casting  a  shadow  over  the  opal  diffuser  of  the  instrument. 

To  test  the  efficacy  of  Equation  24  with  the  Monte  Carlo  simulations,  Ed(Sun)  is  computed 

from 

Ej (Sun)  =  cos  i?o  D0  exp  [  -  (t^  +  tr  4-  toz)/  cos  d0] ,  (25) 

where  Fq  is  the  extraterrestrial  solar  irradiance,  and  r*,  tr,  and  tqx  are,  respectively,  the  contri¬ 
butions  to  the  optical  thickness  of  the  atmosphere  from  aerosol  scattering,  molecular  (Rayleigh) 
scattering,  and  Ozone  absorption.  Ed(Sky)  is  then  determined  by  subtraction  from  the  total  irradi- 
ance  falling  on  the  sea  surface.  Even  though  Equation  25  is  exact,  for  our  purposes  it  underestimates 
Ed(Sun)  because  all  photons  scattered  by  the  aerosol  are  assumed  to  be  uniformly  distributed  over 
the  sky,  whereas  in  reality  a  significant  fraction  of  the  aerosol  scattering  is  through  small  angles  and 
these  scattered  photons  are  still  traveling  in  nearly  the  same  direction  as  the  unscattered  photons. 
To  compensate  for  this  effect,  we  can  obtain  an  upper  limit  on  Ed(Sun)  by  ignoring  the  aerosol 
scattering  entirely,  i.e.,  by  computing  Ed(Sun)  according  to 

Ed(Sun)  =  cos  d0  E0  exp  [  -  (tr  -|-  r0t)l  cos  d0] ,  (26) 

which  clearly  overestimates  Ed(Sun)  since  aerosol  scattering  does  contribute  something  to  Ed(Sky). 
Thus,  for  our  purposes  Equations  25  and  26,  respectively,  provide  lower-  and  upper-bound  estimates 
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of  Ed{ Sun)  and  therefore  of  /.  Comparison  between  Do  computed  from  Equation  24  using  Equation 
25  for  Dd(Sun)  and  the  “exact”  values  (Table  3)  show  that  for  0  <  t?o  <  60°  the  error  is  less  than 
±3%,  and  for  i?o  —  80°  Equation  24  yields  a  value  for  Do  which  is  5-8%  too  low.  The  corresponding 
computations  with  Equation  24  using  Equation  26  for  E<j(Sun)  show  that  for  0  <  t?o  <  60°  the 
error  is  less  than  ^2%,  and  for  t?o  =  80°  the  computed  value  is  0.5-4%  too  high. 

We  can  apply  this  computation  to  the  “shadow”  method  suggested  above  for  estimating  /. 
Assume  that  the  object  used  to  cast  the  shadow  of  the  sun  is  a  circular  disk  of  diameter  somewhat 
larger  than  the  collecting  face  of  the  irradiance  meter.  Then,  if  the  disk  is  relatively  close  to  the 
irradiance  meter,  a  portion  of  the  sky  in  the  vicinity  of  the  sun  is  also  obscured.  This  would 
approximately  correspond  to  estimating  /  using  Equation  26,  i.e.,  photons  scattered  at  small 
angles  from  the  sun  would  be  included  in  Ej(Sun).  Conversely,  if  the  disk  were  at  a  great  distance 
from  the  instrument  only  the  solar  disk  itself  would  be  obscured,  and  photons  scattered  at  small 
angles  from  the  sun  become  part  of  Ej(Sky),  approximately  corresponding  to  using  Equation  25 
to  estimate  /.  Thus  we  conclude  that  the  shadow  method  of  determining  /  should  yield  values  of 
Dq  between  the  estimated  obtained  using  Equations  25  and  26. 

In  the  presence  of  surface  waves,  computation  of  the  correct  value  of  D0  is  facilitated  by  the 
empirical  observation  that  Do  increases  approximately  in  proportion  to  a 2  for  wind  speeds  up  to 
~  20  m / s.  This  is  demonstrated  in  Figure  9  for  an  overcast  sky  and  for  solar  illumination  (no 
atmosphere)  with  i?0  =  60°,  70° ,  and  80° .  The  dots  on  Figure  9  are  the  computed  values  of  D0 
and  the  lines  are  least-squares  fits  to  Do  =  ci  +  C2<r*,  where  Cj  and  cj  are  constants.  The  least- 
squares  lines  allow  the  estimation  of  D0  with  and  error  of  <  2%.  For  i90  <  50°  the  variation  in  D0 
for  0  <  <r  <  0.3  is  less  than  about  2%.  Thus,  for  i?o  £  50°  Do  can  be  computed  by  assuming  that 
the  sea  surface  is  flat,  while  for  larger  values  of  i?  (or  for  an  overcast  sky)  the  flat-surface  values  of 
Do(Sun)  and  Do(Sky)  for  use  in  Equation  23  must  be  increased  in  accordance  with  Figure  9. 

Finally,  in  the  atmospheric  model  used  here  ta  at  550  nm  was  taken  to  be  0.25.  This  is 
very  conservative,  since  it  would  correspond  to  a  coastal  atmosphere  (it  is  typical  of  a  continental 
aerosol)  and  is  a  factor  of  2—3  higher  than  would  be  expected  for  a  “clear”  marine  atmosphere. 
Therefore  an  estimate  of  /  based  on  Equation  26  alone,  i.e.,  without  any  measurements  above  the 
surface,  will  provide  excellent  estimates  of  Dq  in  clear  marine  atmospheres. 
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The  Lambert- Beer  Law  Applied  to  Kd 


Having  established  that  K  and  (K)  can  be  transformed  into  inherent  optical  properties  in  a 
variety  of  realistic  situations,  we  now  turn  to  the  main  question  of  this  paper:  the  extent  to  which 
Kd  satisfies  the  Lambert-Beer  law.  Consider  an  ocean  consisting  of  m  components  one  of  which  is 
pure  sea  water.  Let  c*  be  the  specific  attenuation  coefficient  of  constituent  i.  Then,  c,-  =  c*Cj,  and 
the  total  attenuation  coefficient  can  be  written 

m  m 

c  =  Ec‘  =  Ec*c<’  (2?) 

«=i  »=i 

where  C{  is  the  concentration  of  the  iih  constituent.  These  relationships  comprise  the  Lambert-Beer 
law,  i.e.,  the  individual  attenuation  co  fficients  are  proportional  to  the  individual  concentrations 
and  the  total  attenuation  coefficient  is  a  linear  sum  of  the  individual  or  partied  attenuation  co¬ 
efficient.  [On  the  surface,  Equations  9  emd  10  appear  to  suggest  that  cp  is  not  proportional  to 
the  concentration  of  phytoplankton,  but  rather  on  the  concentration  to  the  0.6  power.  However, 
this  is  an  eirtifact  because  cp  in  the  present  bio-optical  model  includes  not  only  the  contribution  of 
phytoplankton  but  also  the  contribution  from  detrital  material,  the  relative  concentration  of  which 
varies  with  the  concentration  of  phytoplankton  (Hobson,  Menzel  and  Barber  1973).  In  Teality  the 
attenuation  coefficient  of  particles  in  Case  1  waters  should  be  written  cp  =  c*hCph  +  c’dCd,  where 
the  subscripts  “ ph ”  and  “d”  refer  to  phytoplankton  and  detritus,  respectively.]  Since  K/ D0  and 
(K)/ Do  are  inherent  optical  properties,  the  relevant  question  concerns  the  validity  of  the  expres¬ 
sions 

=  ^  Ei 

Do  ~  Do 

and 

(JO  _  V  (JOi 

Do  fri  Do  ' 

For  a  given  observation,  the  Do's  cancel  from  Equations  28  and  29;  however,  we  will  keep  D0  on 
both  sides  of  these  equations  because  Equations  18-21  express  K/Do  and  (K)/Do  as  functions  of 
the  inherent  optical  properties,  and  also  because  later  we  will  consider  combining  measurements 
made  under  differing  environmental  conditions,  i.e.,  measurements  with  different  values  of  Do  as 
would  be  carried  out  in  the  field. 
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Clearly,  if  Equation  21  is  used  for  Kf  Do,  its  linear  dependence  on  the  inherent  optical  prop¬ 
erties  means  that  the  error  in  Equation  28  is  no  more  than  the  error  in  Equation  21,  i.e., 

jy  to  m  to  to  .. 

—  =  1.0395(a  +  6b)  =  1.0395  +  =  £  1  0395(af  +  (&b)f)  =  Y,~nL- 

0  »=1  i=l  »=1  i= 1  0 

Thus,  any  errors  in  Equations  28  and  29  over  and  above  the  error  in  Ki/Do  and  (K)i/Do  will 
result  from  nonlinearities  in  the  dependence  of  these  quantities  on  the  inherent  optical  properties. 
To  understand  the  magnitude  of  these  additional  errors  we  consider  a  hypothetical  model.  Assume 
that  Equation  19,  the  more  nonlinear  of  the  two  relationships,  is  exact  and  that  the  ocean  consists 
only  of  water  and  plankton.  With  uiw  and  wp  given  in  Table  2  and  Fw  =  0.50  and  Fp  =  0.985,  the 
relative  concentration  of  particles,  as  measured  by  ep/c,  is  varied  from  0  to  1.  The  true  value  of 
{ K)/Dq  ,  ( K)t/Dq  ,  is  then  computed  from  Equation  19  using  the  value  of  u0 F  for  the  mixture, 
i.e.,  using 

r,  t^pF pCp  -f-  ujpj FyjC^ 

uqF  =  - - . 

Cp  *T 

The  Lambert-Beer  law  value,  ( K)g/ Dq  ,  is  computed  from 

{K)b  _  (K) u,  (K)p 
Do  Dq  Do 

with  { K)wjD0  and  {K)pj D0  individually  determined  by  Equation  19  using  wwFw  and  u>pFp,  re¬ 
spectively.  Graphically,  in  Figure  7  { K)b/cDo  would  fall  on  a  straight  line  between  the  points  on 
the  least-squares  curve  at  u>o F  =  uiwFw  and  uqF  =  u)pFp,  while  { K)t/cDo  would  be  on  the  curve. 
The  relative  error  in  ( K)b/Do  is  then  computed  by  means  of 

( K)b  -  ( K)t 

(K)t  • 

This  error  is  shown  in  Figure  10  as  a  function  of  the  fraction  of  particles  (cp/c).  It  is  seen  that 
the  maximum  error  at  440  nm  is  ~  3%,  while  the  maximum  error  at  550  nm  is  ~  6%.  Had 
wp(440)  =  0.80  been  used  in  this  example  the  error  at  440  nm  would  have  been  <  2%.  Since 
wp(440)  =  0.86  is  near  the  upper  limit  for  phytoplankton,  e.g  ,  Bricaud,  Morel  and  Preiur  (1983) 
measured  u>p(440)  =  0.88  for  the  coccolithophore  Emiliana  huxleyi  which  is  known  to  be  a  very 
strong  scatterer,  it  is  believed  that  the  error  at  440  nm  will  usually  be  less  than  2%.  The  dotted 
curve  on  Figure  10  corresponds  to  a  mixture  of  water  and  nonabsorbing  particles  (such  as  white 
sand)  at  wavelengths  >  600  nm  for  which  uiw  ~  0.  This  curve  provides  the  maximum  error 
incurred  in  (K)B. 
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Using  Equation  13,  cp/c  can  be  related  to  the  pigment  concentration  for  this  two-component 
example.  At  440  nra,  C  >  0.05  mg/m3  yields  cp/c  >  0.75,  while'at  550  nm  C  >  0.5  mg/m3 
yields  cp/c  >  0.75.  Thus,  for  pigment  concentrations  in  the  range  0.5  <  C  <  10  mg/m3,  Figure 
10  and  the  discussion  above  suggest  that  ( K)b/Do  is  about  1  -  2%  too  low  at  440  nm  and  about 
6%  too  low  at  550  nm.  Furthermore,  considering  the  small  errors  in  Equations  18  and  19,  these 
computations  suggest  that  the  Lambert-Beer  law  for  K*  as  expressed  in  Equations  28  send  29  is 
unlikely  to  be  in  error  by  more  than  5%  at  440  nm  or  10%  at  550  nm. 

It  should  be  noted  that  the  example  above  of  the  error  in  the  Lambert-Beer  law  applied  to 
( K)/Do  concerns  adding  strong  scatterers  (plankton)  to  a  strongly  absorbing  medium  (sea  water). 
If  weak  scatters  (wp  «  0)  are  added  to  sea  water,  the  Lambert-Beer  law  will  be  much  better  satisfied. 
For  example,  if  pure  absorbers  are  added  to  pure  water,  the  maximum  error  in  the  Lambert-Beer 
applied  to  (K)/ D0  will  be  <  0.25%.  The  extension  of  these  computations  to  medie  with  m  >  2  is 
straightforward. 

It  is  important  to  understand  that  the  near-validity  of  the  Lambert-Beer  law  rests  squarely  on 
the  near-linearity  of  the  relationships  shown  in  Figures  6  and  7,  i.e.,  that  the  quantities  involved 
must  be  inherent  optical  properties  is  a  necessary  but  not  sufficient  condition  for  the  validity  of  the 
law.  For  example,  if  all  particles  in  the  ocean  were  sufficiently  small  to  scatter  light  with  the  same 
phase  function  as  pure  sea  water,  the  dependence  of  K/cDo  and  ( K)/cD0  on  l-u/0 F  would  be  given 
by  Figure  8.  In  such  a  case,  if  nonabsorbing  particles  were  mixed  with  strongly  absorbing  water, 
according  to  the  Lambert-Beer  law  the  K's  for  the  resulting  mixture  would  fall  {dong  a  straight 
line  from  1  —  ujqF  =  0.5  to  1  —  u>qF  =  1  (F  —  0.5),  while  the  actual  K' s  would  fall  along  the  curve. 
Clearly,  large  departures  from  Lambert-Beer  law  would  be  seen  for  all  values  of  cp/c  in  such  an 
ocean.  Thus,  the  near-validity  of  the  Lambert-Beer  law  in  the  case  of  a  realistic  ocean  is  seen  to 
result  from  the  interplay  of  three  independent  facts:  (1)  the  dependence  of  the  diffuse  attenuation 
coefficients  on  the  geometric  structure  of  the  light  field  can  be  removed  (division  by  Z?0);  (2)  pure 
sea  water  is  a  much  better  absorber  than  scatterer  at  optical  frequencies  (1  —  Fwuw  >  0.85);  and 
(3)  the  phase  functions  for  particles  suspended  in  the  ocean  differs  significantly  from  that  of  pure 
sea  water  (Figure  1). 

Finally  it  is  of  interest  to  determine  the  accuracy  with  which  one  can  estimate  the  diffuse 
attenuation  coefficient  of  an  ocean  consisting  of  pure  sea  water  alone  through  extrapolation  of  Ad- 
values  measured  in  a  real  ocean  to  the  limit  of  zero  particle  concentration.  As  mentioned  in  the 
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introduction,  this  is  the  scheme  that  Smith  and  Baker  and  others  use  to  estimate  K&  for  pure  water 
(Baker  and  Smith  1982,  Smith  and  Baker  1978a,  Smith  and  Baker  1978b,  Smith  and  Baker  1981). 
For  this  purpose  we  have  computed  (K)  as  a  function  of  c  at  480  nm  by  letting  C  in  Equation  13 
range  from  0  to  4.5  mg/m3.  Figure  11  shows  the  results  for  t?  =  0°,  60°,  and  for  overcast  skies. 
The  lines  on  the  graph  correspond  to  linear  least-squares  fits  to  the  computed  points  with  C  >  0 
(cp  >  0  or  c  >  Cu,),  i.e.,  the  point  on  each  line  corresponding  to  pure  water  was  left  out  of  the  fit. 
The  least-squares  line  was  then  extrapolated  to  cp  =  0  to  determine  (K)  in  the  absence  of  particles. 
This  corresponds  to  extrapolating  C  to  zero  pigment  concentration.  As  seen  from  the  figure,  the 
extrapolated  line  falls  very  close  to  the  computed  values  of  ( K )  for  pure  sea  water.  In  fact,  the 
difference  between  the  computed  rind  extrapolated  values  of  { K)w  are,  respectively,  3.8,  1.9,  and 
1.5%.  =  0°,  60°  ,  and  an  overcast  sky. 

In  this  example,  the  incident  illumination  is  the  same  for  each  value  of  c  along  the  least- 
squares  line.  In  practice  this  would  be  impossible  to  arrange  experimentally.  In  the  field,  each 
data  point  would  likely  correspond  to  a  different  incident  light  field.  However,  we  have  seen  that 
division  ( K )  by  D0  removes  the  effects  of  the  geometric  structure  of  the  light  field.  To  assess  the 
efficacy  of  determining  (K)w  from  extrapolation  to  cp  =  0  in  more  realistic  situations,  for  each 
wavelength  the  above  extrapolation  procedure  was  applied  to  ( K)/D0  obtained  from  all  of  the 
simulations,  i.e.,  all  illumination  conditions  were  treated  equally  and  included  in  the  analysis. 
Figure  12  shows  the  results  of  the  extrapolation  at  480  nm,  and  Table  5  compares  the  extrapolated 
value  of  ( K)w/Do  with  the  true  value  of  (K)w  —  the  value  of  ( K )  for  an  ocean  composed  of  pure 
sea  water  computed  with  the  atmosphere  removed  and  with  the  sun  at  the  zenith.  Table  5  suggests 
that  the  extrapolation  procedure  can  yield  the  true  value  of  {K)w  to  within  about  2%.  However, 
Figure  12  shows  that  large  errors  ir-  [K)w  are  possible  if  it  is  determined  from  a  small  amount  of 
data  for  which  cp  »  cw.  For  example,  if  the  highest  value  of  ( K)/Do  at  c  ss  0.27  m~l  and  the 
lowest  value  at  c  a  0.42  m-1  the  extrapolated  value  of  (K^/Do  would  be  «  0.038  m_1,  an  error  of 
nearly  a  factor  of  2.  Thus,  experimental  determination  of  {K)w  must  be  carried  out  by  excluding 
turbid  waters  from  the  analysis. 
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Concluding  Remarks 


By  simulating  the  transport  of  radiation  in  a  realistic  ocean-atmosphere  system  and  treating 
the  results  as  experimental  data  obtained  under  carefully  controlled  conditions,  it  has  been  shown 
that  K  and  (K)  when  modified  through  division  by  Dq  are,  to  a  high  degree  of  accuracy,  inherent 
optical  properties.  A  simple  scheme  for  estimating  Do  for  individual  experimental  situations  is 
provided.  Furthermore  it  is  shown  that  K/Do  and  (K)/Dq  satisfy  the  Lambert-Beer  law  to  a 
reasonable  degree  of  accuracy,  with  the  largest  error  («  15%  using  {K)]Dq  )  arising  from  mixing 
nonabsorbing  particles,  e.g.,  white  sand,  with  strongly  absorbing  water  (A  >  600  nm);  however, 
leaving  Do  out  of  Equations  1,28,  and  29  will  result  in  extra  variance  in  K&,  K,  and  ( K )  ,  which 
Table  3  suggests  will  be  as  much  as  ±8%  even  for  measurements  restricted  to  0  <  t?o  <  40°. 
In  the  case  of  a  two  component  system  composed  of  pure  sea  water  and  plankton,  the  error  in 
the  application  of  the  Lambert-Beer  law  to  (K)/Do  is  ss  5%  at  440  nm  and  10%  at  550  nm. 
This  accounts  for  the  success  of  Equation  1  for  in-situ  observation  and  analysis  of  phytoplankton 
absorption.  The  near-validity  of  the  Lambert-Beer  law  in  this  situation,  where  there  are  compelling 
reasons  to  believe  that  it  should  fail,  is  traced  to  three  independent  facts:  (1)  the  dependence  of 
the  diffuse  attenuation  coefficients  on  the  geometric  structure  of  the  light  field  can  be  removed;  (2) 
pure  sea  water  is  a  much  better  absorber  than  scatterer  at  optical  frequencies;  and  (3)  the  phase 
functions  for  particles  suspended  in  the  ocean  differs  significantly  from  that  of  pure  sea  water.  If 
any  of  these  facts  were  false  the  Lambert-Beer  law  would  fail.  Finally,  it  is  shown  that  extrapolation 
of  K/ D0  and  ( K)/D0  to  the  limit  c  -»  cw  yields  quantities  which  are  within  2%  of  Kw  and  (K)^  , 
i.e.,  the  value  of  K  and  ( K )  that  would  be  measured  for  an  ocean  consisting  of  pure  sea  water  with 
the  sun  at  the  zenith  and  the  atmosphere  removed. 

The  analysis  of  oceanic  properties  using  is  attractive  because  of  the  relative  simplicity 
of  the  instrumentation  required  for  its  measurement.  The  near- validity  of  the  Lambert-Beer  law, 
for  all  but  the  most  strongly  scattering  of  natural  waters,  allows  the  partial  diffuse  attenuation 
coefficients  (in  the  sense  of  Equations  1,  28,  and  29)  to  be  determined.  Since  the  partial,  as  well 
as  the  total,  functions  ( K  and  ( K )  )  are  in  the  first  approximation  proportional  to  a  +  6b,  for 
those  species  for  which  a  >  6b,  e.g.,  phytoplankton  and  dissolved  organic  material,  measurement 
of  K  or  (K)  provides  a  direct  estimate  of  a  from  Ki/D0  or  ( K)i/D0  .  Thus,  until  the  development 
of  an  in-situ  spectral  absorption  meter,  measurement  of  K j  would  appear  to  be  the  only  available 
in- situ  means  of  estimating  a(A).  Note,  however,  that  in  general  the  medium  will  contain  more 
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than  two  components,  e.g.,  water,  plankton,  and  detritus,  and  the  separation  of  the  components 
can  only  be  carried  out  in  a  statistical  sense. 

Application  of  these  results  to  field  experiments  present  several  difficulties.  The  first,  steins 
from  the  fact  that  K ,  which  satisfies  the  Lambert-Beer  law  better  than  (K)  ,  is  very  difficult  to 
measure  in  practice  due  to  the  strongly  fluctuating  irradiance  at  the  surface  resulting  from  the 
presence  of  surface  capillary  waves,  and  the  difficulty  of  accurately  determining  the  depth  of  the 
instrument  near  the  surface  due  to  the  presence  of  surface  gravity  waves.  Thus,  measurement  of 
{ K )  ,  which  is  significantly  less  influenced  by  the  surface  effects,  is  preferred  from  an  experimental 
point  of  view;  however,  in  the  case  of  oceanic  water,  the  mixed  layer  must  be  sufficiently  deep  so 
that  zlo  =  Tio/c  is  within  the  mixed  layer  and  the  water  may  be  treated  as  homogeneous.  For  the 
limiting  case  of  an  ocean  free  of  particles,  this  would  require  a  mixed  layer  of  «  125,  115,  and  35 
m,  respectively  at  440,  480,  and  550  nm.  A  second  difficulty  concerns  the  determination  of  D 0  in 
the  presence  of  broken  clouds.  In  this  case  Equation  24  will  not  apply  and  the  only  viable  method 
of  is  to  photograph  with  a  fisheye  camera.  Finally,  the  presence  of  whitecaps  on  the  sea  surface 
will  further  modify  the  internal  geometry  of  the  light  filed  and  influence  Do-  Their  effect  cannot 
be  discussed  further  without  knowledge  of  their  optical  properties. 

Unfortunately,  the  results  of  this  paper  are  not  direcfly  applicable  to  measurements  of  Kj  made 
with  moored  instruments  in  a  fixed  configuration  since  K4  is  a  function  of  depth  and  is  dependence 
on  depth  depends  on  c.  [The  dependence  of  K&  on  depth  was  not  relevant  to  the  present  work 
because  measurements  were  considered  to  be  made  at  the  unique  depth,  z  =  0,  or  at  a  given  fraction 
of  surface  irradiance.]  The  question  of  interpretation  of  measurements  from  moorings  will  require 
further  study. 
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Figure  Captions 


Figure  1.  Phase  functions  for  particles  and  water.  To  facilitate  plotting,  phase  functions  “M” 
and  “C”  have  been  multiplied  by  2  and  4,  respectively. 

Figure  2.  Total  phase  function  at  480  nm  as  a  function  of  the  particle  concentration.  Progressing 
from  bottom  to  top  on  the  left  of  the  graph,  Cp/e*,  =  0,  1,  3,  10,  and  100. 

Figure  3.  Computed  dependence  of  K^/c  on  depth  at  440  nm  for  particle  phase  function  “M." 
The  three  cases  are,  from  right  to  left,  for  cf/cw  =  0, 1.4,  and  5.6,  respectively. 

Figure  4.  K/c  as  a  function  of  1  -  u>0.  The  symbol  code  is  provided  on  Figure  1.  Note,  1  -  o»o  = 
a/c. 

Figure  5.  ( K)/c  as  a  function  of  1  —  w<j.  The  symbol  code  is  provided  on  Figure  1.  Note, 
1  —  u0  =  a/c. 

Figure  6.  K/cDq  as  a  function  of  1  —  u>o F.  The  symbol  code  is  provided  on  Figure  1. 

Figure  7.  { K)/cDq  as  a  function  of  1  -  u0F.  The  symbol  code  is  provided  on  Figure  1. 
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Figure  8.  K/cDq  for  a  particle  free  ocean  (Rayleigh  scattering).  Open  circles  are  for  K/cDo, 
while  solid  circles  are  for  (K)/cDo  • 

Figure  9.  D0  as  a  function  of  <r2.  The  curves  from  bottom  to  top  correspond,  respectively,  to  a 
completely  overcast  sky  and  to  solar  illumination  with  =  60°  ,  70° ,  and  80°  . 

Figure  10.  Relative  error  (%)  in  {K)b  as  a  function  of  the  relative  concentration  of  particles. 

Figure  11.  ( K )  at  480  nm  as  a  function  of  c  for  particle  phase  function  “M."  The  lower  and  upper 
lines  are  for  d0  =  0°  and  60°  ,  respectively,  the  center  line  is  for  an  overcast  sky. 

Figure  12.  ( K)/D0  as  a  function  c  at  480  nm.  The  points  are  Monte  Carlo  simulations  for  various 
values  of  t?0,  the  line  is  a  least-squares  fit  to  the  points  with  cp  >  cw. 
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Table  1:  Absorption  and  Scattering  Coefficients  of  Pure  Sea  Water. 


A 

(nm) 

bw 

(m-1) 

aw 

(m-1) 

440 

0.0049 

0.0145 

0.253 

480 

0.0034 

0.0176 

0.162 

550 

0.0019 

0.0638 

0.029 

Table  2:  Model  values  of  u>p 

and  w„. 

A 

Up 

uw 

(nm) 

440 

0.86 

0.253 

480 

0.88 

0.162 

550 

0.93 

0.029 

Table  3:  D0  just  beneath  the  sea  surface. 


t?0 

440  nm 

480  nm 

550  nm 

1/  COSt?0tu 

0° 

1.0343 

1.0267 

1.0185 

1.0000 

20° 

1.0737 

1.0652 

1.0553 

1.0346 

25° 

1.0880 

1.0768 

1.0669 

1.0544 

30° 

1.1050 

1.0995 

1.0933 

1.0787 

o 

O 

1.1580 

1.1540 

1.1487 

1.1415 

o 

O 

<0 

1.2855 

1.2925 

1.2993 

1.3154 

o 

O 

CO 

1.2841 

1.3105 

1.3457 

1.4838 

Diffuse 

1.1969 

1.1969 

1.1969 

Table  4: 
and  •do  - 

Computed  Do  and  diffuse  attenuation  coefficients  at  480  nm 
-  60°  as  a  function  of  the  surface  roughness  parameter  <r. 

<T 

Do 

K/c 

K/cDo 

(K)/c 

{K)/cDo 

0.0 

1.293 

0.2624 

0.2029 

0.2914 

0.2254 

0.1 

1.306 

0.2632 

0.2015 

0.2924 

0.2261 

0.2 

1.333 

0.2733 

0.2050 

0.2954 

0.2285 

0.3 

1.373 

0.2833 

0.2063 

0.2999 

0.2319 

Table  5:  The  true  { K)w  and  the  extrapolated  value  of  (K)w/Do 
in  m_1  for  the  three  wavelengths. 


440  nm  480  nm  550  nm 


(KU 


(K)w/D0 


0.0182 

0.0178 


0.0202 

0.0202 


0.0652 

0.0667 


Figure  1.  Phase  functions  for  pax' 
phase  functions  “M”  and  “C”  have 


((..isis)  q/(©)d  NoixoNnd  asvHd 
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0.4 


Kd/c 

0.6 


0.8 
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Figure  3.  Computed  dependence  of  Kd/c  on  depth  at  440  nm  for  particle 
phase  function  “M.”  The  three  cases  are,  from  right  to  left, 
for  Cp/cu  =  0,  1.4,  and  5.6,  respectively. 


Figure  6.  K, 
The  symbol 


Figure  7.  (K)/cD0  as  a  function  of  1  -  w0F. 
The  symbol  code  is  provided  on  Figure  1. 
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Figure  8.  K/cD0  for  a  pi 
Open  circles  are  for  K/> 


(%)  1{>I>  /  ( X(5I>  -  a(5I» 


Solid:  440  nm 

Dashed:  550  nm 

Dotted:  cop  =  1.00,  cow  =  0.00 


Figure  11.  (K)  at  480  nm  as  a  function  of  c  for  pa 
The  lower  and  upper  lines  are  for  =  0°  and  60° ,  respective!; 


(j.m)  °a/(M> 


Figure  12,  ( K)fD0  as  a  function  c  at  480  nm.  The  poi 
for  various  values  of  tf0,  the  line  is  a  least-squares 
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Abstract 


Through  Monte  Carlo  simulation  the  variation  of  the  diffuse  reflectance  of  natural  waters  with 
sun  tingle  is  found  to  be  dependent  on  the  shape  of  the  volume  scattering  function  (VSF)  of  the 
medium.  It  is  shown  that  single  scattering  theory  can  be  used  to  estimate  the  reflectance  -  sun 
angle  variation  given  the  VSF,  and  conversely,  the  VSF  can  be  retrieved  from  measurements  of  the 
variation  of  the  reflectance  with  sun  angle.  The  complex  variation  of  reflectance  with  the  incident 
illumination  and  surface  roughness  can  be  reduced  to  the  variation  of  a  single  parameter:  the 
downwelling  distribution  function  in  the  absence  of  scattering.  These  observations  are  applicable 
to  all  but  the  most  reflective  of  natural  waters. 


The  irradiance  ratio  or  irradiance  reflectance  at  a  depth  z  is  defined  according  to  R(z)  = 
Eu(z)/Ed(z),  where  Eu  and  Ej  are,  respectively,  the  upwelling  and  downwelling  irradiances  at  z. 
When  this  is  evaluated  just  beneath  the  sea  surface  it  is  referred  to  as  the  diffuse  reflectance  and 
indicated  by  R.  R  is  interesting  because  it  is  relatively  easy  to  measure,  e.g.,  absolute  calibration 
of  the  irradiance  meter  is  not  required,  and  it  is  used  in  the  theory  of  ocean  color  remote  sensing 
(Gordon  and  Morel  1983).  In  a  series  of  Monte  Carlo  simulations  of  the  transport  of  optical 
radiation  in  the  ocean  Gordon,  Brown  and  Jacobs  (1975)  computed  R  for  a  homogeneous  ocean 
as  a  function  of  the  inherent  optical  properties  of  the  water,  the  absorption  coefficient  a,  the 
scattering  coefficient  b  and  the  volume  scattering  function  /3(a),  using  scattering  phase  functions 
[P(a)  =  0(a)/b]  measured  by  Kullenberg  (1968)  in  the  Sargasso  Sea.  In  a  limited  number  of  cases 
fZ(t?o),  the  diffuse  reflectance  as  a  function  of  the  solar  zenith  angle  t?0)  was  also  studied.  It  was 
concluded  that  JZ(i? 0)  was  a  very  weak  function  of  i?o  -  varying  by  less  than  20%  for  0°  <  t?0  <  60°. 
Later,  Kirk  (1984)  presented  a  similar  Monte  Carlo  study  (using  phase  functions  measured  by 
Petzold  (1972)  in  San  Diego  Harbor)  which  showed  a  variation  in  jR(t?0),  over  the  same  range  of 
angles,  of  much  as  50%.  This  difference  is  far  greater  than  what  might  be  expected  due  to  differences 
in  the  computational  procedure  and  thus  requires  explanation.  If  it  can  be  verified  that  neither 
computation  is  in  error,  differences  in  the  results  can  only  lie  in  the  fact  that  different  scattering 
phase  functions  were  used  in  the  computations.  In  fact,  Jerlov  (1976,  page  149)  states  that  the 
variation  of  R  with  t?o  “is  a  consequence  of  the  shape  of  the  scattering  function”,  but  provides  no 
quantitative  demonstration  of  the  claim.  Also,  even  for  i90  =  0,  Plass,  Kattawar  and  Humphreys 
(1981)  have  already  shown  that  R  depends  on  the  shape  of  the  scattering  phase  function,  contrary 
to  the  conclusion  of  Gordon,  Brown  and  Jacobs  (1975),  Kirk  (1984),  and  Morel  and  Prieur  (1977) 
that  R  depends  on  the  phase  function  principally  through  the  backscattering  coefficient  bf,  given  by 

bf,  =  2nb  /  P{ a)  sin  a  da ; 

Jn/ 2 

however,  they  employed  phase  functions  that  differed  considerably  from  those  observed  in  natural 
waters  in  order  to  demonstrate  the  dependence. 
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To  investigate  quantitatively  the  influence  of  the  scattering  phase  function  on  the  diffuse  re¬ 
flectance,  I  have  recalculated  R  using  standard  Monte  Carlo  techniques  and  also  using  a  new 
backward  Monte  Carlo  code  developed  for  another  purpose  (Gordon  1985).  The  results  of  the  new 
computations  are  in  complete  agreement  with  the  old,  i.e.,  when  the  Petzold  (1972)  phase  function 
is  used  the  results  from  either  code  agree  well  with  Kirk  (1984)  and  when  the  Kullenberg  (1968) 
phase  function  is  used  the  results  agree  with  Gordon,  Brown  and  Jacobs  (1975).  This  suggests  the 
computations  of  both  Kirk  (1984)  and  Gordon,  Brown  and  Jacobs  (1975)  were  correct  and  that 
the  difference  is  in  fact  due  to  the  specific  phase  functions  used  in  the  two  studies.  It  also  suggests 
that  the  shape  of  the  variation  of  R  with  i?o  can  provide  some  information  about  /3(a) .  In  what 
follows,  it  is  shown  that,  given  /3(a),  the  single  scattering  approximation  can  be  used  to  specify  the 
variation  of  R  with  i?o,  and  conversely,  given  R(d o),  it  is  possible  to  invert  the  process  and  estimate 
(3(a). 

Let  the  ocean  be  illuminated  by  the  direct  solar  beam.  Then,  if  is  the  solar  zenith  angle 
observed  beneath  the  sea  surface,  i.e.,  msin  t?o»  =  sint?o,  where  m  is  the  refractive  index  of  water, 
in  the  single  scattering  approximation  the  diffuse  reflectance  is  given  by 

(1) 

where  cos  a  =  ~npw  +  y/(l-  p2)(l  -  /i*,)cos  <p,  and  =  cosdow  This  expression  is  valid  for 
6  <  a,  and  provides  a  direct  link  between  R(0O)  and  0(a)  in  this  limit.  Gordon  (1973)  has  shown 
that,  for  a  medium  that  scatters  strongly  in  the  forward  direction,  the  validity  of  Equation  1  can  be 
extended  somewhat  by  replacing  c  by  ct  +  Physically  this  corresponds  to  replacing  c  by  the  value 
of  the  downwelling  irradiance  attenuation  coefficient  (Ka)  that  would  be  measured  just  beneath 
the  surface  of  the  medium  (with  the  sun  at  the  zenith),  and  it  is  called  the  quasi-single  scattering 
approximation  (QSA).  Typically  for  the  ocean  b t,  <C  6,  and  the  QSA  holds  for  much  larger  values 
of  b/a  (or  equivalently  u>0  =  b/(a  +  b))  than  the  single  scattering  approximation.  With  the  sun 
at  the  zenith,  Gordon  (1973)  found  the  the  QSA  reproduced  Monte  Carlo  computations  of  R  to 
within  0.5%  for  u>o  <  0.6  and  12%  for  uq  <  0.85.  Thus,  one  can  expect  Equation  1  to  predict 
the  variation  of  R  with  sun  angle,  even  for  rather  large  values  of  u/0.  The  efficacy  of  Equation 
1  in  this  respect  is  examined  by  comparison  with  exact  (Monte  Carlo)  computations  carried  out 
using  two  scattering  phase  functions:  “KC"  ,  the  Kullenberg  (1968)  phase  function  measured  at 
460  nm  in  the  Sargasso  Sea  and  used  by  Gordon,  Brown  and  Jacobs  (1975);  and  “T”,  the  mean 
of  the  three  particle  phase  functions  measured  in  turbid  water  at  530  nm  by  Petzold  (1972)  and 
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used  in  the  Kirk  (1984)  computations.  These  are  shown  in  Figure  1  along  with  phase  function  for 
molecular  (Rayleigh)  scattering  of  pure  water.  The  resulting  comparison  between  Equation  1  and 
the  Monte  Carlo  (M.C.)  computations  is  presented  in  Figure  2,  where  it  is  seen  that  the  analytical 
computation  of  R(tf0)/R( 0)  agrees  with  the  Monte  Carlo  computations  with  a  maximum  error  of 
<  5%  for  wo  =  0.8  and  <  10%  for  u>o  =  0.9.  The  general  increase  in  R  with  i?o  is  due  to  the  fact 
that  the  minimum  scattering  angle  for  photons  to  be  redirected  to  the  surface  increases  with  i?o 
according  to  a  -  t/2-i?ou>-  Since  P(a)  increases  rapidly  with  decreasing  a  for  a  <  w/2  (Figure  1), 
portions  of  the  phase  function  which  are  larger  than  those  for  t/2  <  a  <  it  increasingly  contribute 
to  R  as  i?0  increases.  This  is  particularly  evident  for  phase  function  T  for  which  there  is  a  strong 
contrast  in  the  values  of  P(a)  above  and  below  a  =  t/2.  In  the  case  of  scattering  by  pure  sea 
water  (not  shown),  i.e.,  Rayleigh  scattering,  the  analytical  computation  has  a  maximum  error  of 
<  0.5%  for  the  largest  value  of  w0  encountered  in  the  ocean  (~  0.3  near  400  nm). 

It  is  seen  that  for  w0  =  0.8,  there  is  a  clear  trend  for  the  multiple  scattering  computations 
of  J2(i?o)/R(0)  to  be  above  the  single  scattering  results.  This  is  due  to  the  fact  that  at  each 
scattering  event  after  the  first,  some  photons  can  scatter  through  progressively  smaller  angles  and 
still  be  redirected  toward  the  surface.  Since  P(a)  increases  with  decreasing  a  for  a  <  t/2,  multiple 
scattering  should  increase  iZ( i?o)/ P(0)  above  the  single  scattering  result  when  photons  can  scatter 
a  few  times  before  being  absorbed.  Since  phase  function  T  is  larger  at  very  small  scattering  angles 
than  KC  (25%  of  the  scattering  events  for  T  have  a  <  1°,  compared  to  5%  for  KC),  the  effect  is 
larger  for  KC.  For  larger  values  of  wo,  this  effect  disappears  because  photons  scatter  many  times 
before  reaching  the  surface  and  the  information  concerning  the  direction  at  which  the  photon 
entered  the  water  becomes  lost.  Thus,  for  small  w0  the  exact  value  of  P(i?o)/P(0)  should  be  close 
to  that  given  by  Equation  1;  however,  as  wq  increases,  J2(i?o)/I2(0)  will  initially  increase  above  the 
single  scattering  value  and  then  eventually  decrease  below  it  as  wo  — ►  1. 

The  near-agreement  between  P(t?o)/P(0)  and  that  predicted  by  single  scattering  (or  more  cor¬ 
rectly,  QSA)  indicates  that  measurements  of  /3(a)  over  the  appropriate  range  of  angles,  41°  <  a  < 
180  ,  can  be  used  to  predict  the  dependence  of  R  on  i?oi  however,  measurement  of  /3(a)  over  this 
entire  range  is  not  necessary.  Gordon  (1976)  showed  that  64  could  be  accurately  determined  with¬ 
out  knowing  the  full  scattering  function.  This  is  accomplished  by  fitting  measurements  of  / 3  at  only 
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three  angles,  a  =  45°,  90°  ,  and  135° ,  to  an  analytic  equation  first  used  by  Beardsley  and  Zaneveld 
(1969): 


/3(a)  = 


/3(9Q°) 


'  '  ~  (1  -  e/  cosa)4(l  +  efccosa)4  ’ 

where  e f  and  eb  are  adjustable  parameters.  The  fits  of  phase  functions  T  and  KC  to  Equation  2 
are  shown  as  the  solid  lines  in  Figure  3.  The  fit  is  excellent  for  scattering  between  about  40°  and 
160° ,  correctly  reproducing  the  significant  variation  around  90° .  However,  Equation  2  is  a  very 
poor  approximation  at  scattering  angles  less  than  25°  -30°  and  is  in  error  by  a  factor  of  103  or 
more  near  0° .  For  Rayleigh  scattering  (not  shown),  Equation  2  provides  an  excellent  fit  for  all 
scattering  angles.  Computation  of  P(do)/P(0)  using  Equation  1  (single  scattering)  and  the  fits  of 
the  phase  function  to  Equation  2  agree  with  those  using  Equation  1  and  the  actual  phase  functions 
with  an  error  of  less  than  2%.  On  this  basis  I  conclude  that  measurement  of  /3(45°),  /?( 90°),  and 
/3(135°)  are  sufficient  to  describe  the  variation  of  R  with  do.  It  is,  of  course,  of  interest  to  know 
if  the  process  above  can  be  inverted,  i.e.,  given  P(do)/P(0)  is  it  possible  to  estimate  /9(a)?  To 
examine  this  question,  I  have  assumed  that  R( do)  i*  measured  (in  this  case  simulated  by  Monte 
Carlo)  at  10°  increments  from  0  to  89° .  This  “data”  is  then  fit  to  Equation  1,  with  /3(a)  given 
by  Equation  2,  using  a  nonlinear  least-squares  technique  to  determine  the  unknown  parameters  ej 
and  eb.  The  resulting  derived  phase  functions  for  u>o  =  0.8  (dashed  curves)  and  0.9  (chain  dashed 
curves)  are  presented  in  Figure  3.  In  reality,  the  above  procedure  can  only  provide  P(a)/P(90°); 
however,  in  Figure  3  we  have  normalized  P(a)  to  obtain  the  correct  value  of  bb,  i.e.,  so  that  the  true 
P(a)  and  the  inverted  P(a)  would  yield  the  same  value  of  P(0).  Clearly,  in  this  restricted  case,  i.e., 
measurement  of  fZ(d0)  for  the  full  domain  of  do,  the  general  shape  of  / 3(a)  for  60°  —  70°  <  a  <  180° 
can  be  retrieved  from  measurements  of  P(d o). 


The  conclusions  in  this  note  are  based  on  simulations  of  an  idealized  ocean,  i.e.,  a  flat,  homo¬ 
geneous  ocean  with  wo  <  0.9  in  the  absence  of  the  atmosphere.  How  applicable  are  they  to  a  real 
ocean?  From  the  analysis  of  Gordon  (1987)  one  expects  uo  <  0.9  in  natural  waters  except  in  intense 
plankton  blooms  (in  the  green)  or  regions  with  a  high  concentration  of  nonabsorbing  suspended 
particles,  e.g.,  white  sand,  in  the  water.  The  only  effect  of  the  atmosphere  is  to  add  a  quasi-difluse 
component  (skylight)  to  the  irradiance  incident  on  the  sea  surface.  If  dEj  is  the  irradiance  incident 
on  the  surface  with  zenith  angles  between  d  and  d  +  <fd,  the  reflectance  is  easily  computed  to  be 

R  !o/2  P(t?)T(d)(dPrf/dd)dd 
f0*/2  T(4)(dEd/d4)d4  ’ 
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where  T(d)  is  the  Fresnel  transmittance  (air  towater)  of  the  surface  for  an  incident  angle  of 
t?.  For  a  totally  diffuse  incident  irradiance  (incident  radiance  independent  of  viewing  direction) 
d Ed /d-d  <x  sin  2d  and  R  normalized  to  R( t?o  =  0)  is  approximately  that  given  by  Equation  1  for 
d0  =  40°  in  agreement  with  Monte  Carlo  simulations.  Similarly,  the  effect  of  replacing  the  fiat 
surface  with  a  rough  surface  is  to  render  the  incident  light  field  beneath  the  surface  more  diffuse. 
However,  in  contrast  with  a  flat  ocean,  i?oto  can  exceed  the  critical  angle,  sin_1(l/m),  when  the 
surface  is  rough.  For  large  do  large  values  of  i?o«>  are  possible,  and  these  photons  can  can  be 
redirected  toward  the  surface  through  highly  probable  small  angle  scattering.  Thus,  in  the  absence 
of  skylight  i?(i?0)  will  in  general  be  larger  in  an  ocean  ruffled  by  the  wind  than  for  a  flat  ocean. 
Monte  Carlo  simulations  using  a  surface  described  by  the  Cox  and  Munk  (1954)  slope  distribution 
characteristic  of  a  7.2  m/s  wind  speed  with  phase  function  T  and  u?o  =  0.8  showed  slight  (<  2%) 
increases  in  R(do)/R(0)  over  the  smooth  ocean  case  for  t?o  =  40°  and  60°  and  increases  of  7  and  18% 
for  t?o  =  70°  and  80°  ,  respectively.  Thus  the  effect  of  surface  waves  on  the  validity  of  R(do)f  R(0) 
computed  from  Equation  1  is  small  for  <  60  -  70°. 

It  would  be  useful  if  the  effect  on  R  of  variations  in  t?o,  the  amount  of  skylight  in  the  incident 
irradiance,  and  the  surface  roughness  could  be  explained  by  the  variation  of  a  single  parameter.  The 
key  to  finding  such  a  parameter  lies  in  a  recent  set  of  Monte  Carlo  simulations  (Jerome,  Bukata  and 
Buxton  1988)  showing  that  iZ(i?o)/iJ(0)  =  1/  cos  t?ow  Although  the  present  computations  only  have 
the  property  that  R(d0)/R(Q)  «  kf  cosi?0u>,  where  &  is  a  constant,  the  value  of  ik  is  approximately 
unity  when  the  same  phase  function  used  by  Jerome,  Bukata  and  Burton  (1988)  is  employed.  The 
values  of  k  for  phase  functions  KC  and  T  are  approximately  0.85  and  1.15,  respectively.  These 
observations  are  more  interesting  when  it  is  realized  that  if  the  illumination  incident  on  a  flat 
ocean  is  in  the  form  of  a  parallel  beam  from  the  sun,  the  downwelling  distribution  function  —  the 
downwelling  scalar  irradiance  divided  by  the  downwelling  irradiance  (Preisendorfer  1961)  —  just 
beneath  the  surface  in  the  absence  of  scattering,  Do,  is  exactly  l/cost?o».  In  cases  with  parallel 
beam  illumination  of  a  flat  ocean,  as  well  as  situations  with  more  complex  illumination,  Do  has  been 
used  to  remove  the  geometrical  properties  of  the  incident  light  field  from  the  downwelling  irradiance 
attenuation  coefficient  (Gordon  1989,  Gordon,  Brown  and  Jacobs  1975),  in  effect  normalizing  the 
coefficient  to  that  which  would  be  measured  in  the  absence  of  the  atmosphere  with  the  sun  at  the 
zenith.  Thus,  it  is  natural  to  ask  if  Do  might  be  used  to  simplify  the  analysis  of  R(do)/R(0)  in 
more  complex  situations.  Clearly,  the  simplest  procedure  for  trying  to  extend  the  analysis  to  more 
complex  situations  is  to  plot  R(d0)/R(0)  as  a  function  of  Do,  i.e.,  to  replace  the  1/cos t?0»  used 
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by  Jerome,  Bukata  and  Burton  (1988)  by  its  more  general  equivalent,  Do-  Thus  far,  in  addition 
to  an  incident  beam  of  parallel  irradiance  on  a  flat  surface,  we  have  discussed  two  other  incident 
distributions  (beneath  the  surface):  a  distribution  resulting  from  wind  induced  surface  roughness; 
and  a  distribution  resulting  from  uniform  radiance  incident  above  a  flat  surface.  In  the  case  the 
wind-ruffled  surface  described  earlier,  a  separate  computation  was  carried  out  to  determine  Do  and 
the  resulting  values  of  5(Z?o)/-S(l)  are  plotted  against  Do  in  Figure  4  (symbols  with  +'s)  along  with 
the  results  taken  from  Figure  2.  Note  that  in  each  case,  the  simulated  values  of  R(Do)/R(l)  are 
linearly  related  to  D0,  and  that  the  rough  ocean  cases  fall  with  excellent  accuracy  along  the  same 
lines  as  their  flat  ocean  counterparts.  Simulations  carried  out  with  totally  diffuse  incident  irradiance 
falling  on  a  flat  ocean  (symbols  with  x's  in  Figure  4)  show  that  the  same  linear  relationship  is 
satisfied  in  this  case  as  well. 

These  results  show  that  for  a  given  u>0  and  scattering  phase  function,  R(D0)  =  kD0R(l), 
i.e.,  the  variation  of  R  with  the  incident  illumination  and  the  surface  roughness  can  be  completely 
explained  through  their  effect  on  Do-  [An  accurate  scheme  for  estimating  Do  from  simple  irradiance 
and  wind  speed  measurements  is  provided  in  Gordon  (1989).]  The  parameter  k  depends  mostly 
on  the  scattering  phase  function,  and  for  u>o  <  0.9  it  can  be  computed  using  the  single  scattering 
approximation.  In  sum,  the  main  conclusion  of  this  note  —  that  single  scattering  can  be  used  to 
characterize  the  variation  of  R  with  the  incident  radiance  distribution  —  should  be  applicable  to  a 
real  ocean  if  highly  reflective  waters  ( u>o  >  0.9)  are  avoided. 
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Figure  Captions 


Figure  1.  Phase  functions  used  in  the  present  study. 

Figure  2.  R  as  a  function  of  t?o-  Dots  are  for  u>o  =  0.8,  open  circles  for  u>o  =  0.9.  The  solid 
curves  are  the  result  of  the  single  scattering  approximation  (Eq.  1).  Upper  curve  is 
for  phase  function  “T”,  lower  for  “KC.” 

Figure  3.  Phase  functions  “KC”  (upper)  and  “T”  (lower).  The  dots  are  the  true  values.  The 
solid  lines  are  fits  with  Eq.  2  using  P(ot)  evaluated  at  a  =  45°,  90° ,  and  135° .  The 
dashed  and  chain  dashed  curves  are  the  result  of  inverting  Eq.  1  for  cu0  =  0.8  and  0.9, 
respectively. 

Figure  4.  R  as  a  function  of  Do.  Dots  are  for  u>o  —  0.8,  open  circles  for  w0  =  0.9.  Symbols  with 
+'s  are  for  a  rough  surface.  Symbols  witE  x's  are  for  a  diffuse  incident  irradiance. 
The  solid  curves  are  the  result  of  the  single  scattering  approximation  (Eq.  1).  Upper 
curve  is  for  phase  function  “T,”  lower  for  “KC.” 
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Figure  1.  Phase  functions  used  in  the  present  study. 
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In  a  recent  paper1  we  investigated  the  effects  of  multiple 
scattering  on  the  Coastal  Zone  Color  Scanner  (CZCS)2-4 
atmospheric  correction.  As  a  measure  of  the  efficacy  of  the 
correction  algorithm,  we  defined  the  residual  error  at  443  nm 
in  the  correction  algorithm  AL„(443),  when  the  water  leaving 
radiance  was  known  in  the  other  three  bands  (520,  550,  and 
670  nm),  according  to 

A£c(443)  =  L'( 443)  -  £,(443)  -  r(443)£J443) 

FJ 443) 

-<(443.670)  f£.(670)-r 

where  L;(X),  L,(X),  and  Lw(\)  are,  respectively,  the  total 
radiance  observed  by  the  sensor  at  the  wavelength  X,  the 
radiance  due  to  Rayleigh  scattering  in  the  atmosphere,  and 
the  radiance  exiting  the  sea  surface  from  beneath;  t (X)  is  the 
diffuse  transmittance  of  the  atmosphere,  and  Fo(X)  is  the 
extraterrestrial  solar  irradiance.  <(443,670)  is  the  so-called 
atmospheric  correction  parameter  that  is  determined  from 
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Fig.  1.  Error  in  the  recovered  aerosol  radiance  at  443  nm  (in  digital  counts)  across  the  CZCS  scan  for  Orbit  2217  as  a  function  of  the 
stratospheric  aerosol  concentration  (r,)  and  phase  function:  (a)  t„(670)  *0.2  with  the  Haze  L  phase  function  in  the  stratosphere;  (b)  to(G70)  » 
0.4  and  Haze  L  in  the  stratosphere;  (c)  r„(670)  *■  0.2  and  Haze  M  in  the  stratosphere;  (d)  t„(670)  -  0.4  and  Haze  M  in  the  stratosphere. 


spectral  variation  of  the  radiance  resulting  from  aerosol  scat¬ 
tering  in  the  atmosphere  and  is  approximately  the  ratio  of 
the  aerosol  optical  thickness  at  443  to  that  at  670  nm.  For  a 
perfect  atmospheric  correction,  ALa(443)  =»  0;  for  an  accept¬ 
able  atmospheric  correction  AL<,(443)  5  2  DC,  where  1  DC 
(digital  counts)  is  the  radiance  corresponding  to  the  digitiza¬ 
tion  interval  (one  part  in  256)  of  L,(X).  At  443  nm,  1  DC  ®= 
0.045  mW/cm2pm  sr.  We  computed  AL„(443)  for  the  CZCS 
orbit  as  a  function  of  the  scan  angle  and  the  aerosol  optical 
thickness  in  the  red  and  found  that  the  simple  expediency  of 
using  the  standard  algorithm,5  which  was  justified  on  the 
basis  of  single  scattering,  computing  the  molecular  scatter¬ 
ing  contribution  by  correctly  accounting  for  multiple  scatter¬ 
ing,  and  reducing  *he  satellite-derived  value  of  {(443,670)  by 
5%,  reduced  AZ„(443)  to  below  the  1-2  DC  level  for  most 
atmospheres  and  viewing  situations.  In  those  computations 
the  aerosol  was  largely  confined  to  the  troposphere.  Howev¬ 
er,  the  eruption  of  the  volcano  El  Chichon  during  28  Mar.-4 
Apr.  1982  was  known  to  have  introduced  a  significant  quanti¬ 
ty  of  aerosol  into  the  stratosphere.6  Since  radiative  transfer  , 
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in  the  atmosphere  is  not  insensitive  to  the  vertical  distribu¬ 
tion  of  the  components,  and  since  CZCS  viewed  the  ocean 
through  the  El  Chichon  aerosol  layer  for  a  significant  portion 
of  its  life  (1978-1986),  we  believed  it  would  be  useful  to 
extend  the  computations  to  include  the  effect  of  an  El  Chi- 
chon-Iike  aerosol  layer  on  CZCS  atmospheric  correction. 
We  present  the  results  of  these  simulations  in  this  Letter. 

The  tropospheric  model  follows  that  used  previously,1  i.e., 
a  stratified  atmosphere  with  both  the  Rayleigh  and  aerosol 
scattering,  each  decreasing  exponentially  with  altitude,  but 
with  different  scale  heights.  Approximately  90%  of  the  aero¬ 
sol  is  confined  to  a  layer  of  2  km  in  thicknes*  near  the  sea 
surface,  with  ~80%  of  the  Rayleigh  scattering  molecules 
above  the  aerosol.  To  simulate  the  effects  of  the  El  Chichon 
stratospheric  aerosol,  a  homogeneous  layer  of  arbitrary  opti¬ 
cal  thickness  r,  and  phase  function  is  added  at  the  top  of  the 
atmosphere.  The  tropospheric  and  stratospheric  aerosols 
both  scatter  according  to  two-term  Henyey-Greenstein7 
(TTHG)  phase  functions.  Two  aerosol  models  are  used  in 
the  computations.  The  first  approximates  marine  aerosol 


phase  functions  given  by  Quenzel  and  Kastner.8  The  TTIIG 
parameters  for  this  phase  function,  which  we  shall  refer  *•.  as 
Haze  XI,  are  a  =  0.983.  =  0.82.  and g->  =  —0.55.  1  he  second 

phase  function  is  an  approximation  developed  by  Kattawar7 
ft  r  the  Deirmendjian  Haze  L  distribution9  with  a  refractive 
index  of  1.55,  used  to  represent  continental-type  aerosols. 
Its  parameters  are  a  =  0.9618, g\  =  0.7130,  and#2  =  —0.7598. 
In  our  simulations,  the  marine  aerosol  model  is  always  used 
for  the  tropospheric  phase  function.  Interestingly,  this 
phase  function  does  not  differ  significantly  in  the  backward 

i.^.. .sphere  from  that  proposed  for  the  El  Chichon  strato- 
spheric  aerosol  by  King8  based  on  extinction  measure¬ 
ments.1^  11  Thus,  in  one  set  of  simulations,  the  stratospheric 
aerosol  phase  function  is  also  approximated  by  the  marine 
aerosol  model.  In  a  second  set  of  simulations  the  Haze  L 
phase  function  is  used  for  the  stratospheric  aerosol.  Haze  L 
differs  considerably  front  King’s  model,  scattering  approxi¬ 
mately  an  order  of  magnitude  more  at  scattering  angles  O 
near  180°.  Haze  L  probably  represents  the  maximum  scat¬ 
tering  near  0  =>  180°  for  realistic  aerosol  models.  It  was 
included  to  provide  contrast  to  King’s  model.  The  tropo¬ 
spheric  aerosol  optical  thickness  is  assumed  to  vary  with 
wavelength  according  to  -a(A)  =  (670/A)nra(670),  where  n  =  0 
or  1,  while  the  stratospheric  aerosol  optical  thickness  is  as¬ 
sumed  to  be  independent  of  wavelength.10  11  The  small 
effects  of  ozone  are  ignored. 

Computations  for  the  four  orbital  scenarios  studied  earli¬ 
er1  were  carried  out  with  the  stratospheric  aerosol  layer  in 
place.  A  sample  of  the  computations  is  given  in  Fig.  1,  which 
provides  ALU(443)  at  nine  positions  across  the  CZCS  scan12 
as  a  function  of  ra(670)  and  ts  for  orbit  2217.  The  computed 
AL  .1443)  values  are  for  r„(670)  =  0.20  and  0.40,  n  =  1,  and  r, 
=  0,  0.10,  and  0.20.  Using  the  standard  definition  of  visible 
range,  i.e..  3.912 /6(550>,  where  6(550)  is  the  total  scattering 
coefficient  at  the  surface  at  550  nm,  a  ra(670)  of  0.20  with  n  = 
1  corresponds  to  a  horizontal  visible  range  at  the  surface  of 
-“15.4  km  while  0.40  corresponds  to  7.9  km.  The  Haze  L 
phase  function  has  been  used  for  the  stratospheric  aerosol  in 
Figs.  1(a)  and  (b),  while  Figs.  1(c)  and  (d)  are  for  the  Haze  M 
phase  function.  The  symbol  •  at  a  computed  point  indi¬ 
cates  that  the  total  radiance  at  the  CZCS  was  greater  than 
the  saturation  radiance  (256  DC)  at  670  nm.  As  the  aerosol 
concentration  increases,  the  sensor  tends  to  saturate  because 
of  the  contribution  by  Rayieigh  scattering  at  large  scan  an¬ 
gles.  At  these  points  no  atmospheric  correction  is  possible  in 
practice  because  the  total  radiance  at  the  sensor  at  670  nm  is 
unknown. 

Figures  1  show  the  effect  of  the  vertical  distribution  of  the 
aerosol  on  AL0(443).  Comparing  the  r,  =  0.20  curve  in  (c) 
[ r,, (670)  =  0.20)  and  r,  =  0  in  pane!  (d)  (ra(670)  =  0.40J 
reveals  an  ~1.5-DC  difference  in  AL„(443);  however,  these 
curves  correspond  to  the  same  total  aerosol  optical  thickness 
(ru(670)  +  r,  =  0.40),  and  both  aerosols  are  characterized  by 
the  same  phase  function  (Haze  M).  The  only  difference 
between  the  two  curves  is  that  in  (c)  half  of  the  aerosol  is  near 
the  bottom  of  the  atmosphere  and  half  of  the  aerosol  is  at  the 
top  „f  the  atmosphere,  while  in  fd)  all  the  aerosol  is  near  the 
bottom  of  the  atmosphere. 

The  immediate  and  most  important  conclusion  to  be 
drawn  from  Fig.  1  is  that,' except  near  the  edges  of  the  scan,  in 
those  situations  where  the  sensor  does  not  saturate  at  670  nm 
the  error  in  the  atmospheric  correction  algorithm  in  the 
presence  of  the  stratospheric  aerosol  is  <1-2  DC.  This 
observation  is  also  applicable  to  the  other  viewing  scenarios 
examined  earlier,1  with  the  exception  of  Orbit  3226  with 
Haze  L  in  the  stratosphere.  The  strong  hackscattering  asso- 
i  iaU-d  with  Haze  U  comes  into  play  in  this  orbit  where  almost 
diret  t  hackscattering  is  observed  near  the  center  of  the  scan 


(O  -  174. h”).  This  causes  the  sensor  to  saturate  even  near 
the  scan  center  with  rs  as  low  as  0.10  when  r„(670)  =  0.40and 
at  several  points  along  the  scan  when  r a(670)  =  0.20.  When 
the  sensor  does  not  saturate  near  the  center  of  the  scan  the 
error  is  <3  DC  for  tq(670)  =  0.20  and  4  DC  for  t„(670)  =  0.40. 
When  the  Haze  M  phase  function  is  used  in  the  stratosphere 
in  Orbit  3226,  a  pattern  similar  to  that  in  Figs.  1(c)  and  (d)  is 
found;  however,  for  ro(670)  =  0.40  and  ts  =  0.20  the  error  is 
slightly  larger  than  2  DC.  It  is  interesting  to  note  that  the 
addition  of  the  stratospheric  aerosol  actually  improves  the 
correction  in  some  cases.  For  example,  Figure  11(c)  in  Ref.  1 
shows  that  near  the  scan  center  of  Orbit  2381  the  error  is 
approximately  -1.5  DC  in  the  absence  of  a  stratospheric 
aerosol,  while  our  calculations  (not  presented)  reveal  that 
with  the  addition  of  a  stratospheric  aerosol  with  r3  =  0  2  this 
error  is  reduced  to  about  -0.25  DC. 

Finally,  we  note  that  the  presence  of  the  aerosol  layer  will 
decrease  the  scene  contrast  over  what  would  be  observed  if  it 
were  all  in  the  troposphere,  because  adjacency  effects,11  i.e., 
the  atmospheric  scattering  of  photons  leaving  the  surface 
from  pixels  other  than  the  one  under  examination,  are  a 
strong  function  of  the  vertical  distribution  of  the  aerosol. 
However,  because  of  the  small  reflectance  of  the  ocean  and 
the  generally  low  contrast  of  ocean  scenes,  this  effect  is 
believed  to  be  unimportant  for  CZCS.14 

In  summary,  the  computations  presented  here  suggest 
that  the  addition  of  an  El  Chichon-like  aerosol  layer  in  the 
stratosphere  has  very  little  effect  on  the  basic  CZCS  atmo¬ 
spheric  correction  algorithm.  The  main  influence  of  the 
additional  stratospheric  aerosol  is  to  increase  the  total  radi¬ 
ance  exiting  the  atmosphere  and  thereby  increase  the  proba¬ 
bility  that  the  sensor  will  saturate.  In  the  absence  of  satura¬ 
tion  the  correction  algorithm  should  perform  as  well  as  in  the 
absence  of  the  stratospheric  layer. 

This  work  received  support  from  the  National  Aeronautics 
and  Space  Administration  under  grant  NAGW-273  and  con¬ 
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For  measurement  of  aerosols  over  ihe ocean.  the  total  radiance  I.,  backscattered  from  the  top  of  a  stratified  at¬ 
mosphere  which  contains  both  stratospheric  and  tropospheric  aerosols  of  various  types  has  been  computed. 
A  similar  computation  is  carried  out  for  an  aerosol-free  atmosphere  yielding  the  Rayleigh  scattered  radiance 
I  he  difference  L:  —  l.r  is  shown  to  be  linearly  related  to  the  radiance  L m,  which  the  aerosol  would  prudiu  e 
in  the  single  scattering  approximation.  This  greatly  simplifies  the  application  of  aerosol  models  to  aerosol 
analysis  by  satellite  since  adding  to,  or  in  some  way  changing,  the  aerosol  model  requires  no  additional 
multiple  scattering  computations.  In  fact,  the  only  multiple  computations  required  for  aerosol  analysis  are 
those  lor  determining  L:,  which  can  be  performed  once  and  for  all.  The  computations  are  explicitly  applied  to 
Hand  1  of  the  I'ZCS.  which,  because  of  its  high  radiometric  sensitivity  and  excellent  calibration,  is  ideal  tor 
studying  aerosols  over  the  ocean.  .Specifically,  the  constant  A  in  the  relationship  L0>  =  A  "  —  /..!  is  given  as 

a  function  ot  position  along  the  scan  for  four  typical  orbital-solar  position  scenarios.  The  computations  show 
that  can  be  retrieved  from  L,  -  L.  with  an  average  error  of  no  more  than  5-7'<  except  at  the  very  edges  of 
the  scan. 


I.  Introduction 

There  is  considerable  interest  in  th*>  global  distribu¬ 
tion  of  aerosols  because  of  their  role  in  climate  and 
biogeochemical  cycling.1  Because  our  knowledge  of 
aerosol  concentrations  in  remote  areas,  particularly 
over  the  oceans,  is  very  limited  due  to  the  paucity  of 
measurements,  there  has  been  considerable  effort  di¬ 
rected  toward  estimating  aerosol  concentration  and 
other  properties  using  earth-orbiting  satellites.  "9 

As  a  result  of  the  potential  success  promised  by  the 
earlier  investigations,  NASA  included  aerosol  radiance 
as  one  of  four  standard  output  products  from  the  Nim- 
bus-T  Coastal  Zone  Color  Scanner  (CZCS).U)  The 
CZCS  is  a  scanning  radiometer  which  views  the  ocean 
in  six  coregistered  spectral  bands,  five  in  the  visible 
and  near  IK  (443,  520,  550,  670,  and  750  nm,  labeled 
Bands  1,  2,  3,  4,  and  5,  respectively)  and  a  thermal  IR 
band  ( 10.5-12.5  pm,  Band  6).  The  sensor  has  an  active 
scan  of  78°  centered  on  nadir  and  a  field  of  view  of 
0.0495°,  which  from  a  nominal  height  of  955  km  pro¬ 
duces  a  ground  resolution  of  825  m  at  nadir.  The 
satellite  is  in  a  sun-synchronous  orbit  with  ascending 
node  near  local  noon.  The  sensor  is  equipped  with 
provision  for  tilting  the  scan  plane  ±20°  from  nadir  in 
2°  increments  along  the  satellite  track  to  minimize  the 
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influence  of  direct  sun  glint.  Th  ?  purpose  of  the  CZCS 
experiment  was  to  provide  estimates  of  the  near  sur¬ 
face  concentration  of  phytoplankton  pigments  (the 
sum  of  the  concentrations  of  chlorophyll  a  plus  phaeo- 
phytin  a)  and  total  seston  by  measuring  the  spectral 
radiance  backscattered  out  of  the  ocean."  The  radi¬ 
ance  backscattered  from  the  atmosphere  and  sea  sur¬ 
face  (specular  reflection)  is  typically  at  least  an  order 
of  magnitude  larger  than  the  desired  radiance  scat¬ 
tered  out  of  the  water.  Therefore,  the  CZCS  has  very 
highradiometricsensitivity.  Overmuch  of  the  world's 
ocean,  the  radiance  backscattered  out  of  the  water  in 
Band  4  is  below  the  sensitivity  limit  of  the  sensor  and 
can,  therefore,  be  taken  to  be  zero.  Thus,  this  combi¬ 
nation — a  very  sensitive  and  well-calibrated  scanning 
radiometer  viewing  a  target  (the  ocean)  which  is  essen¬ 
tially  black — makes  the  CZCS  an  ideal  tool  for  oceanic 
aerosol  studies.  Also,  we  shall  see  that  studying  aero¬ 
sols  with  CZCS  imagery  is  relatively  inexpensive  from 
a  computational  standpoint  since  the  required  quanti¬ 
ty  is  already  generated  at  the  first  step  in  the  atmo¬ 
spheric  correction  procedure.1213 

At  present  the  scheme  for  analysis  of  satellite  imag¬ 
ery  for  aerosols,  as  applied  to  NOAA’s  Advanced  Very 
High  Resolution  Radiometer  (AVHRR)  by  Griggs,4 
involves  the  computation  of  the  expected  radiance  at 
the  sensor  as  a  function  of  the  aerosol  concentration  for 
various  viewing  angles  and  solar  zenith  and  azimuth 
angles.  The  computations  are  then  placed  in  a  look¬ 
up  table,  and  aerosol  concentration  (optical  thickness) 
is  determined  from  the  radiance  measurements  and 
sun-viewing  geometry.  This  requires  a  model  of  the 
optical  properties  of  the  aerosol.  The  model  was  cho¬ 
sen  to  reproduce  previous  measurements  made  with 
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Landsat  2  over  'he  ocean  at  San  Diego.  C A.'-’  In  some 
cases  it  may  be  desirable  to  use  a  different  model  of  the 
aerosol,  e.g.,  the  aerosol  may  be  known  to  be  principal¬ 
ly  composed  of  Saharan  dust.7  or  the  spectral  variation 
of  backscatlering  may  suggest  a  particular  si/e  fre¬ 
quency  distribution. u  This,  however,  requires  an  ex¬ 
tensive  set  of  computation-  for  each  new  aerosol  mod¬ 
el.  In  this  paper  we  present  a  scheme  for  studying 
aerosols  with  CZCS  which  circumvents  having  to  do 
any  radiative  transfer  computations  involving  the 
properties  of  the  aerosol. 

II.  Computation  of  the  Aerosol  Radiance 

The  atmospheric  correction  algorithm  for  CZCS  im¬ 
agery  has  been  described  in  considerable  detail  previ¬ 
ously1’ and  is  not  discussed  here.  Limiting  our  in¬ 
terest  to  CZCS  Band  4  at  670  nm  and  noting  that  the 
water  is  essentially  black  in  this  hand, 1,1  i.e.,  the  water¬ 
leaving  radiance  represents  <1  digital  count  from  the 
sensor,  the  single  scattering  approximation  allows  on 
to  partition  the  total  radiance  at  the  sensor  L,  into  a 
component  due  to  Rayleigh  (r|  scattering  Lrs  and  a 
component  due  to  aerosol  u  scattering 

in 

and  /. ...  are  given  by 

/ . . ,  =  ^xr >  4.t  rusfl.  (21 

where 

p.'-v'j  =  j  +  [>*«">  +  /*<" 

eostf  t  -  ±  cu$4V,  —  sinb,.  co$(o  —  4>,  >■ 

tl,  and  do.  are.  respectively,  the  zenith  and  azimuth 
angles  of  a  vector  from  the  point  on  the  sea  surface 
under  examination  (pixel)  to  the  sun,  and  likewise  0 
and  o  are  the  zenith  and  azimuth  angles  of  a  vector 
from  the  pixel  to  the  sensor.  /.(«)  is  the  Fresnel  reflec¬ 
tance  of  the  interface  for  an  incident  angle  9,  PtU>)  is 
the  scattering  phase  function  of  component  Ad*  =  r  or 
a),  w,  is  the  single  scattering  albedo  of  .v(u.>  =  1 ),  and  tx 
is  the  optical  thickness  of  x(rr  =  0.044).  F„  is  the 
instantaneous  extraterrestrial  solar  irradiance  F<>  re¬ 
duced  by  two  trips  through  the  ozone  layer, 

F,  «  f exp[— r0,(l/ci.s«  +  l/fo.-e,,)). 

where  rU2  is  the  ozone  optical  thickness.  The  term 
involving  8-  in  Eq.  (2)  provides  the  contribution  due  to 
photons  which  are  hackscattered  from  the  atmosphere 
without  interacting  with  the  sea  surface.  The  term 
involving  °  +  accounts  for  those  photons  which  are  scat¬ 
tered  in  the  atmosphere  toward  the  sea  surface  (sky 
radiance)  and  then  specularly  reflected  from  the  sur¬ 
face  into  the  field  of  view  of  the  sensor  termj  as 
well  as  photons  which  are  first  specularly  reflected 
from  the  sea  surface  and  then  scattered  by  the  atmo¬ 
sphere  into  the  field  of  view  of  the  sensor  (/>(0„)  term). 
Clearly,  [,a,  is  directly  related  to  the  optical  properties 
of  the  aerosol  through  uapa(»,8h)  and  the  aerosol  col¬ 
umn  concentration  through  Given  the  optical 
properties  of  aerosol,  its  mass  concentration  can  he 


estimated  from  Las.  Unfortunately,  because  of  multi¬ 
ple  scattering,  Eq.  ( 1)  is  incorrect  and  must  he  replaced 
by1317 

(31 

where  Lr  and  La  are  the  multiple  scattered  counter¬ 
parts  of  Lrs  and  Las ,  respectively,  and  CU  I'  accounts  for 
the  interaction  between  Rayleigh  and  aerosol  scatter¬ 
ing.  CRfl  can  be  either  positive  or  negative  and  in  fact 
can  even  sign  along  a  CZCS  scan  line.13  There  is  no 
way  to  retrieve  La  from  the  other  terms  in  Eq.  (6) 
because  CHP  is  unknown.  Furthermore,  even  if  L„ 
could  be  determined,  unlike  Lai,  it  is  not  simply  related 
to  rQ,  u;0,  and  pa(8,00).  Since  Las  has  a  simple  interpre¬ 
tation  it  seems  desirable  to  be  able  to  extract  it  from  L, 
or  at  least  to  relate  it  to  quantities  extracted  from 
In  this  paper  we  use  radiative  transfer  theory  to  com¬ 
pute  L,  and  Lr  (Lr  =  L,  when  ra  =  0)  including  all  orders 
of  multiple  scattering  in  the  scalar  approximation,  i.e., 
polarization  is  ignored.  Then  the  difference  L,  -  Lr  is 
formed  and  compared  with  Lai.  We  will  show  that  La, 
can  be  retrieved  from  Lt  —  Lr  across  the  CZCS  scan  in  a 
manner  that  is  relatively  insensitive  to  the  ( in  general 
unknown)  scattering  phase  function  of  the  aerosols 
and  thus  derive  a  quantity  Las  which  can  he  directly 
related  to  aerosol  models  without  having  to  pass 
through  the  transfer  equation. 

III.  Atmospheric  Model  and  Computations 

The  atmospheric  model  follows  that  used  in  Cordon 
and  Castano.13  Briefly,  it  consists  of  a  stratified  atmo¬ 
sphere  w;th  the  Rayleigh  scattering  coefficient  b,  de¬ 
creasing  exponentially  with  altitude  z  according  to  a 
scale  height  Hr  and  the  aerosol  coefficient  bn  also  de¬ 
creasing  exponentially  with  altitude  with  a  scale  height 
Ha.  i-e., 

b,  =  b°  expl-Z''//,),  <41 

where  x  and  r  or  a,  and  Hr  and  Ha  are  taken  to  he  9.2 
and  1  km,  respectively.18  Assuming  no  absorption,  the 
Rayleigh  and  aerosol  optical  thicknesses  are  given  by 

r  x  =  b*H,. 

Thus  ~90%  of  the  aerosol  is  confined  to  a  layer  of  2-km 
thickness  near  the  sea  surface  with  of  the  Ray¬ 
leigh  scattering  molecules  above  the  aerosol.  This 
should  be  typical  of  oceanic  aerosols  generated  at  the 
sea  surface  and  confined  for  the  most  part  to  the  ma¬ 
rine  boundary  layer19  and  naturally  generated  conti¬ 
nental  aerosols.20  To  account  for  the  possibility  of 
stratospheric  aerosol  layers,  e.g.,  created  by  volcanic 
activity,  a  homogeneous  aerosol  layer  of  arbitrary  opti¬ 
cal  thickness  r,  and  phase  function  P,{0)  is  placed  at 
the  top  of  the  atmosphere.  The  tropospheric  and 
stratospheric  aerosols  scatter  according  to  a  two-term 
Henvey-Greenstein  phase  function 

/\<«>  =  ,1  -Ml  -  <.)/(«,«,).  15) 

where 
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/> "JC 1  = - -  - —7  • 

(1  +  H'  -  -A’  i  orf'l '  ■ 

and  x  =  a  or  x.  Three  aerosol  models  are  used  in  the 
computations.  The  f  irst  approximates  marine  aerosol 
phase  functions  given  by  Quenzel  and  Kastner.21  It 
has  been  used  by  Sturm--  tor  atmospheric  cot  rection  of 
CZt  'S  imagery  and  for  ocean  color  sensitivity  analyses. 
The  parameters  in  Eqs.  lb)  and  (6)  are  a  =  0.983, #1  = 
0.82,  and  -  — O.bb.  The  second  is  a  two-term  Hen- 
yt  y-Greenstein  app.  jximation  developed  by 
Kattawar2'  to  the  Deirmendjian  Haze  1,  distribution-1 
with  a  refractive  index  of  l.bb  and  is  used  to  repre  ,ent 
continental-!.'  pe  aerosols.  Its  parameters  are  «•  = 
0.9018,  AO  =  0.7310,  and  go  =  —0.7598.  The  third  is  a 
two  term  Henyev-Greenstein  approximation  to  tne 
Deirmendjian  Haze  C  distribution  with  a  slope  pa¬ 
rameters  v  -  3.5  and  a  refractive  index  of  1.50.  Its  pa¬ 
rameters  are  a  =  0.9618,  £1  =  0.7130,  and  go  =  —0.50. 
These  three  phase  function?  are  shown  in  Fig.  1,  where 
the  Haze  L  case  has  the  largest  hackscattering  (0  = 
180°),  Haze  C  the  intermediate  backscattering,  and 
the  marine  aerosol  model  has  the  lowest  backscatter¬ 
ing.  In  most  cases  the  aerosol  is  assumed  to  be  nonab¬ 
sorbing;  however,  a  few  cases  with  =  0.9  have  also 
been  examined. 

The  ocean  is  assumed  to  be  flat  and  totally  ebsorb- 
ing:  i.e.,  all  photons  tha  penetrate  the  surface  are  ab¬ 
sorbed.  The  transport  equation  is  solved  by  the  meth¬ 
od  of  successive  orders  of  scatte.ing  (see,  e.g.,  Ref.  25), 
and  the  results  are  then  transformed  into  CZCS  scan 
coordinates.  The  computations  for  a  given  atmo¬ 
sphere.  i.e..  a  given  set  Pu,  Ps,  t0.  r.„  and  rr,  are  first 
cai  ried  out  to  compute  L,.  Then  r„  and  r,  are  set  to 
zero,  and  the  trai  sport  equation  is  again  solved  to 
determine  Lr.  In  this  manner  L,  -  Lr  can  be  found  for 
a  given  atmosphere-orbital  geometry  scenario. 

A  total  of  four  orbital  scenarios  is  examined,  two 
with  no  sensor  tilt  and  two  with  a  forward  tilt  of  20°. 
These  scenes  are  listed  by  orbit  number  in  Table  I. 
Orbit  130  corresponds  to  fall  in  the  Gulf  of  Mexico, 
while  Orbit  2381  corresponds  to  spring-fall  viewing 
near  the  Arctic  Circle  (Iceland).  These  are  typical  of 
orbits  for  which  the  sensor  is  operated  in  the  nontilted 
mode.  Such  operation  is  marginal  for  Orbit  130,  which 
shows  some  evidence  of  sun  glitter  near  the  center  of 
the  scan  in  the  southern  portion  of  the  Gulf.  The 
other  two  orbits  2217  and  3226  correspond  to  viewing 
in  the  spring  and  summer  in  the  Middle  \tlantic  Bight, 
where  the  low  solar  zenith  angle  necessitates  tilting  the 
scan  forward  to  avoid  the  sun  glitter.  The  geometry  of 
Orbit  3226  was  similar  to  that  which  existed  for  the 
■  ZCS  validation  study  presented  by  Gordon  et  al.lb 
The  variation  in  the  local  solar  zenith  angle  across  the 
scan  (column  6  in  Table  I)  is  surprisingly  small,  consid¬ 
ering  that  the  scan  swath  is  ~  1600  km  on  the  ground  in 
the  untilted  mode  and  2400  km  with  a  tilt  of  20°.  In 
contrast  for  the  tilted  scan,  ft 0  at  the  western  edge  is 
significantly  larger  than  for  the  rest  of  the  scan.  The 
nominal  value  of  (/0  used  in  the  computations  for  each 
orbit  is  given  in  column  7.  The  local  viewing  angle  ft 


Fig.  1.  Aerosol  phase  functions  used  in  the  study:  ▼  ,  Haze  L;  A. 
Haze  C;  •,  marine  aerosol  model. 


varies  between  0°  at  the  center  of  the  scan  and  .66° 
at  the  edges  of  the  scan  in  the  untilted  cases  and 
between  23.16  and  61.01“  in  the  eases  with  a  20°  tilt. 
Thus  in  the  t’lted  mode  the  sensor  views  the  ocean 
through  an  air  mass  that  changes  by  nearly  a  factor  of  2 
over  the  scan. 

IV.  Results 

First,  consider  cases  for  which  there  is  no  aerosol  in 
the  stratosphere,  i.e.,  Tj=  9.  We  computed  L,  -  Lr  and 
Las  for  the  three  aerosol  phase  functions  in  Fig.  1,  the 
orbital  geometries  listed  in  Table  I,  and  four  aerosol 
optical  thicknesses,  r„  =  0.05, 0.20, 0.40,  and  0.60.  A  ra 
of  0.20  is  somewhat  more  turbid  than  the  V23  (23-km 
visibility)  model  used  by  Viollier  et  al.,11 26  while  ra  = 
0.60  is  approximately  the  same  value  used  in  their  V5 
(5-km  visibility)  model.  Since  our  i  laze  C  phase  func¬ 
tion  corresponds  to  that  used  by  Deschamps  et  al.,11 
our  results  for  that  phase  function  should  span  the 
range  of  their  computation,  i.e.,  surface  visibilities 
down  to  ~5  km. 

In  our  preliminary  analysis  of  the  resulting  compu¬ 
tations  we  discovered  that  orbits  with  similar  sensor 
tilts  produced  similar  relationships  between  L;  -  Lr 
and  La$.  Thus  in  our  presentation  we  combine  Orbits 
130  and  2381  with  a  tilt  of  zero  and  Orbits  2217  and 
3226  with  a  tilt  of  20°.  Figures  2  and  3  provide  Las  as  a 
function  of  Lt  —  Lr  at  the  center  of  the  CZCS  scan  for 
the  untilted  and  tilted  orbitals,  respectively.  The 
units  for  the  radiances  are  CZCS  Gain  1  digital  counts 


Table  I.  Scan**  Examined  In  the  Present  Study* 


Orbit 

Tilt 

Lat.6 

Long.* 

DavA'r 

80 

80  (Nominal) 

130 

0 

28 

-86.1 

306/1978 

42.4-43.9 

43.0 

2381 

0 

66 

-22.9 

104/1979 

58.0-58.8 

58.3 

2217 

20 

38 

-62.2 

92/1979 

37.7-39.3r 

39.0 

3226 

20 

38 

-64.4 

165/1979 

19.1-22. 7“* 

22.3 

“  All  angles  are  in  degrees. 

b  Lat.  and  long,  refer  to  6ie  suborbital  latitude  and  longitude  at 
nadir. 

'  40.9  at  the  western  edge  of  the  scan. 
d  25.4  at  the  western  edge  of  the  scan. 
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Fig  2  La,  as  a  function  of  Lt  -  L,  at  t fie  center  of  the  scan  for  the 
nuntilted  orbits  with  the  aerosol  restricted  to  the  troposphere.  •, 
▼  .  and  A  refer,  respectively,  to  the  marine  aerosol,  Haze  L,  and  Ha2e 

C 


Fig  3.  La,  as  a  function  of  L,  -  L,  at  the  center  of  the  scan  for  the 
tilted  orbits  (tilt  =  20°)  with  the  aerosol  restricted  to  the  tropo¬ 
sphere.  •  ,  ▼  ,  and  a  refer,  respectively,  to  the  marine  aerosol,  Haze 
L,  and  Haze  C. 


(DC), where  1  DC  =  0.01103 mW/cm2^msr.  Recalling 
that  the  output  signal  of  the  CZCS  is  8-bit  digitized 
aboard  the  spacecraft,  the  maximum  sensor  output  is 
255  DC.  For  the  geometries  shown,  Lr  **  50  -*  75  DC, 
so  L(  -  Lr  must  be  less  than  ~175-200  DC  to  prevent 
sensor  saturation.  Computations  for  which  Lt  >  255 
DC  have  been  oipitted  from  the  analysis,  since,  al¬ 
though  correct,  they  do  not  correspond  to  observable 
situations  for  the  present  CZCS.  The  most  striking 
feature  of  Figs.  2  and  3  is  the  strong  linearity  between 
L,  -  Lr  and  La,  regardless  of  the  scattering  phase 
function  chosen  to  represent  the  aerosols  or  the  solar 
zenith  angle  associated  with  the  viewing  situation.  In 
Figs.  4  and  5  a  similar  presentation  is  made  for  a  scan 
mirror  position  so  that  the  sensor  is  viewing  a  ground 
position  246  pixels  from  the  eastern  edge  of  the  scan.27 
Note  that  the  linearity  of  the  Lr  -  Lt  vs  Las  relationship 
still  appears  valid.  However,  the  slope  is  larger  than  at 
the  scan  center.  Also,  for  the  tilted  orbits  (Figs.  3  and 


<-o.  (DC) 

Fig.  4.  La,  as  a  function  of  L,  -  L,  at  a  position  located  246  pixels 
from  the  eastern  edge  of  the  scan  ta  scan  mirror  rotation  angle  of 
30°).  The  sensor  is  in  the  untilted  mode,  and  the  aerosol  is  restrict¬ 
ed  to  the  troposphere.  •,  ▼.and  A  refer,  respectively,  to  the  marine 
aerosol,  Haze  L,  and  Haze  C. 


L„,  (DC) 

Fig.  5.  La,  as  a  function  of  L,  —  L,  at  a  position  located  246  pixels 
from  the  eastern  edge  of  the  scan  (a  scan  mirror  rotation  angle  of 
30°).  The  sensor  is  in  the  tilted  mode  (tilt  =  20°),  and  the  aerosol  is 
restricted  to  the  troposphere.  #,  ▼  ,  and  a  refer,  respectively,  to  the 
marine  aerosol,  Haze  L,  and  Haze  C. 


5)  the  Lt  —  Lr  values  appear  to  be  smaller  near  the 
eastern  edge  than  at  the  center  of  the  scan.  This 
illusion  is  caused  by  the  fact  that  three  points  with  Las 
>100  DC  in  Fig.  3  have  caused  the  sensor  to  saturate  at 
the  edge  and  are,  therefore,  not  plotted  on  the  figure. 
This  saturation  is  due  to  an  increase  in  Lr  from  ~75  DC 
at  the  center  to  ~120  DC  at  the  position  in  question 
near  the  eastern  edge  of  the  scan.  For  the  untilted 
orbits  (Figs.  2  and  4)  Lr  is  a  much  weaker  function  of 
position  along  the  scan,  and  the  above  effect  does  not 
occur  with  the  values  of  ra  used  here  (ra  <  0  60).  For 
the  tilted  orbits  the  Haze  L  phase  function  ▼  with  its 
strong  scattering  near  180°  saturates  the  instrument 
with  Ta  as  low  as  0.40. 

One  normally  expects  multiple  scattering  to  increase 
La  over  the  corresponding  Las ;  i.e.,  if  multiple  scatter- 
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ing  computations  were  carried  out  for  the  aerosols 
alone  (r,  =  0),  we  would  expect  La  >  Las.  This  is 
clearly  the  case  in  Figs.  2,  4.  and  5;  however,  in  Fig.  3, 
the  scan  center  for  the  tilted  orbits,  Lt  —  Lr~  Las.  This 
is  the  result  of  the  CRP  term  in  Eq.  (3).  Gordon  and 
Castano13  have  computed  CHP  at  443  nm  (Band  1) 
across  the  CZCS  scan  for  Orbit  2217  for  all  three  aero¬ 
sol  phase  functions  and  t„  =  0.05  and  0.20.  They  find 
that  CRP  >  0  at  all  scan  angles  for  the  marine  aerosol 
model,  CRP  is  slightly  negative  near  the  scan  center 
and  positive  elsewhere  for  the  Haze  C  model,  and  CR'P 
is  large  and  negative  near  the  scan  center  and  slightly 
positive  near  the  scan  edges  for  the  Haze  L  model. 
CR  P  for  Band  4  will  in  general  be  less  than  for  Band  1. 
However,  they  are  qualitatively  similar  so  we  can  ex¬ 
pect  the  same,  albeit  weaker,  behavior  in  Band  4. 
Equation  (3)  shows  that  when  CR  P  is  negative,  it  has 
the  tendency  to  cancel  the  increase  in  La  over  Las 
resulting  from  multiple  scattering.  Thus  in  Fig.  3  the 
marine  aerosol  •  and  the  Haze  C  A  points  tend  to  fall 
above  the  line  Lt  —  Lr  =  Las,  while  the  Haze  L  points  T 
tend  to  fall  below  the  line.  Since  CK  P  is  positive  for  all 
three  phase  functions  when  the  mirror  rotation  angle  is 
30°  east,  this  anomalous  behavior  is  not  seen  in  Fig.  5. 
It  is  interesting  to  note  that  Orbits  130  (tilt  =  0°)  and 
2217  (tilt  =  20°)  have  nearly  the  same  solar  zenith 
angle  (43  and  39°,  respectively)  across  the  scan,  but  the 
dependence  of  Las  on  Lt  —  Lr  is  significantly  different. 
This  is  also  a  manifestation  of  the  CRP  effect.  Thus 
the  La,  vs  Lt  —  Lr  relationship  is  seen  to  depend  on  the 
sensor  tilt  as  well  as  the  solar  zenith  angle  and  the 
position  in  the  scan. 

T o  see  if  the  linearity  of  the  Lt  -  Lr  vs  Las  is  affected 
by  the  addition  of  stratospheric  aerosols,  e.g.,  from  a 
volcano  such  as  El  Chichon,  several  sets  of  computa¬ 
tions  have  been  carried  out  for  cases  with  rs  ^  0. 
Specifically,  the  marine  aerosol  model  was  used  for  the 
tropospheric  phase  function  and  two  values  of  r0  were 
used — 0.20  and  0.40.  Two  models  were  used  for  the 
stratospheric  aerosol:  the  marine  aerosol  phase  func¬ 
tion  and  the  Haze  L  phase  function.  The  values  used 
for  t,  were  0.10  and  0.20.  This  yields  twenty  computa¬ 
tions  of  Lt  -  Lr  for  each  orbit  or  a  total  of  forty  compu¬ 
tations  for  each  sensor  tilt.  The  total  aerosol  optical 
thickness  of  the  atmosphere  is  now  ra  +  rs.  Las  is 
computed  by  summing  the  stratospheric  and  tropo¬ 
spheric  contributions  in  Eq.  (2),  i.e., 

7-0.  -  k'aP,<Mo)  +  .  (7) 

^  4jt  cos# 

As  before  both  aerosols  are  taken  to  be  nonabsorb¬ 
ing.  The  results  of  these  computations  are  presented 
in  Figs.  6-9,  which  are  directly  comparable  with  Figs. 
2-5  for  t,  -  0.  In  Figs.  6-9  the  symbol  ♦  refers  to  cases 
where  r,  ^  0;  however,  for  clarity  the  two  stratospheric 
aerosol  phase  functions  have  not  been  differentiated 
on  the  figures.  Also,  to  try  to  examine  the  effect  of 
particle  absorption  on  the  Lt  -  Lr  vs  Las  relationships, 
three  computations  for  a  mildly  absorbing  (u>a  =  0.90) 
marine  aerosol  in  the  troposphere  (r,  =  0)  were  carried 
out  for  Orbit  2217  (tilt  =  20°).  These  are  indicated  by 
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Fig.  6.  La,  as  a  function  of  L,  -  L,  at  the  center  of  the  scan  for  the 
nontilted  orbits.  Aerosols  are  located  in  both  the  troposphere  and 
the  stratosphere.  •,  ▼.  x^d  ▲  refer,  respectively,  to  the  marine 
aerosol,  Haze  L,  and  Haze  C,  restricted  to  the  troposphere.  ♦  refers 
to  cases  with  aerosols  in  both  the  troposphere  and  the  stratosphere, 
i.e.,  r,  ^  0. 


L„.  (DC) 

Fig.  7.  La,  as  a  function  of  L,  -  Lr  at  the  center  of  the  scan  for  the 
tilted  orbits  (tilt  =  20°).  Aerosols  are  located  in  both  the  tropo¬ 
sphere  and  stratosphere.  •,  ▼,  and  A  refer,  respectively,  to  the 
marine  aerosol,  Haze  L,  and  Haze  C,  restricted  to  the  troposphere. 
♦  refers  to  cases  with  aerosols  in  both  the  troposphere  and  strato¬ 
sphere,  i.e.,  r,  x  0,  and  +  refers  to  cases  for  which  the  tropospheric 
aerosol  has  t'o  =  0.9. 

the  symbol  +  on  the  appropriate  figures.  Comparison 
of  the  corresponding  figures  from  the  two  sets  (2-5  and 
6-9)  show  that  the  Lt  —  Lr  vs  Las  relationship  is  nearly 
unaltered  by  the  presence  of  stratospheric  or  mildly 
absorbing  tropospheric  aerosols. 

We  have  generated  Lt  -  Lr  and  Las  for  the  tropo¬ 
spheric-stratospheric  aerosol  combinations  in  Figs.  6- 
9  at  nine  equally  spaced  points  along  the  CZCS  scan. 
These  correspond  to  mirror  rotation  angles  starting 
40°  west  of  the  subsatellite  track  (-40°)  and  continu¬ 
ing  every  10°  to  40°  east  of  the  subsatellite  track 
(+40°).  At  each  scan  angle  we  fit  La,  by  least  squares 
to  the  linear  relationship 

ALm  =*  (L,  -  Lr )  (8) 

and  determine  A.  The  corresponding  values  of  A,  the 


u 

Q 


S 


T  1 

Ob»3  1  30  and  2381 
Mirror  Rotation  30*  [ail 


% . 


/ 


T 


o'S  .  i  ....  I 

0  50  100 


ISO  200  250 


La.  (DC) 

Fig.  8.  L„,  as  a  function  of  I.,  -  Lr  at  a  position  located  246  pixels 
from  the  eastern  edge  of  the  scan  (a  scan  mirror  rotation  angle  of 
80°).  The  sensor  is  in  the  untilted  mode.  Aerosols  are  located  in 
both  the  troposphere  and  stratosphere.  •,  ▼,  and  A  refer,  respec¬ 
tively,  to  the  marine  aerosol,  Haze  L,  and  Haze  C  restricted  to  the 
troposphere.  ♦  refers  to  cases  with  aerosols  in  both  the  troposphere 
and  stratosphere,  i.e.,  t,  x  0. 


Fig.  q.  La,  as  a  function  of  Lt  —  L,  at  a  position  located  246  pixels 
from  the  eastern  edge  of  the  scan  (a  scan  mirror  rotation  angle  of 
30°).  The  sensor  is  in  the  tilted  mode  (tilt  -  20°).  Aerosols  are 
located  in  both  the  troposphere  and  stratosphere.  •,  ▼,  and  A 
refer,  respectively,  to  the  marine  aerosol,  Haze  L,  and  Haze  C  re¬ 
stricted  to  the  troposphere.  ♦  refers  to  cases  with  aerosols  in  both 
the  troposphere  and  stratosphere,  i.e.,  r,  ^  0,  and  +  refers  to  cases 
for  which  the  tropospheric  aerosol  has  «o  =  0.9. 


255  DC)  and  do  not  represent  a  truly  improved  linear 
relationship.)  Thus,  in  Eq.  (2)  the  combination  (o>a- 
TaPa  +  wsT,ps)  can  be  estimated  along  the  CZCS  scan 
line  from  the  Rayleigh  removed  radiance  L,  -  Lr ■  L,  ~ 
Lr  is  computed  in  the  first  step  of  the  atmospheric 
correction  procedure  and  is  one  of  the  four  standard 
derived  products  produced  from  CZCS  imagery  by 
NASA.10  The  combination  (u>a7apa  +  u),rspj  is  gener¬ 
ally  all  that  can  be  determined  directly  from  CZCS. 
To  proceed  further  other  information  or  an  aerosol 
model  is  required.28  For  example,  one  may  know  from 
lidar  observations  that  r,  ~  0.  Then,  assuming  the 
tropospheric  aerosol  is  nonabsorbing,  a  model  of  the 
aerosol  type,  e  g.,  marine,  Haze  L,  Haze  C,  yields  the 
optical  thickness  ra  at  each  pixel.  It  is  important  to 
note  that  the  various  aerosol  models  can  be  applied 
directly  to  a  derived  quantity  (a >„Tapa  +  w,rsp,)  without 
the  necessity  for  further  radiative  transfer  computa¬ 
tions.  Thus  the  analysis  of  CZCS  imagery  for  aerosol 
concentration  etc.  can  be  carried  out  completely  in  the 
single  scattering  approximation,  even  in  an  atmo¬ 
sphere  with  a  total  optical  thickness  as  high  as  0.644 
(the  largest  value  used  in  our  computations).  This 
considerably  simplifies  the  analysis. 

In  search  of  a  simple  straightforward  way  of  extend¬ 
ing  these  results  to  other  sensor  tilts  and  sun  angles,  we 
tried  to  relate  A  to  the  phase  angle  of  the  observation, 
i.e.,  6 _  in  Eq.  (2).  The  results  showed  that  a  coarse, 
marginally  useful,  relationship  exists  between  A  and 
0~:  A  -  1  +  O.OO7330-,  where  0-  is  in  degrees.  This 
relationship  predicted  values  of  A  for  an  intermediate 
tilt  (12°)  with  a  maximum  error  of  ~7%;  however,  for  a 
tilt  angle  of  20°  and  a  fixed  phase  angle  the  error  it 
induced  in  A  could  be  as  much  as  -7  and  +10%  on  the 
eastern  and  western  sides  of  the  scan,  respectively. 

V.  Summary 

For  measurement  of  aerosls  over  the  ocean,  we  com¬ 
puted  the  total  radiance  Lt  backscattered  from  the  top 
of  a  stratified  atmosphere  which  contains  both  strato¬ 
spheric  and  tropospheric  aerosols  of  various  types.  A 
similar  computation  is  carried  out  for  an  aerosol-free, 
i.e.,  purely  Rayleigh  scattering,  atmosphere  yielding 
the  Rayleigh  radiance  Lr.  The  difference  L,  —  Lr  is 
then  shown  to  be  linearly  related  to  the  radiance  Las, 


number  of  data  points  used  in  the  determination  of  A, 
and  the  average  %  error  in  the  computed  Las  (given  Lt  - 
Lr),  A L„„  defined  By 


AL 


as 


f.„(true)  -  Lai(fit) 
La,(true) 


X  100%, 


are  provided  in  Table  II.  It  is  seen  that  Eq.  (8)  can 
provide  La,  from  Lt  -  Lr  to  within  ~5-8%  except  at  the 
very  edges  of  the  scan.  [Note  that  in  the  case  of  the 
tilted  orbits  the  apparent  improvement  in  the  fit  at 
angles  of  ±40°  over  that  for  ±30°  is  due  to  the  fact  that 
considerably  fewer  points  have  been  used  in  the  fit, 
because  points  were  dropped  due  to  saturation  (L(  > 


Table  It.  Unaar  Regrttalon  Eetlmataa  of  tha  Coafflclant  A  In  Eq.  (8)  aa  a 
Function  of  tha  MhtOf  Rotation  Anglo  a.  A Lu  la  tha  Avarage  Error  In 
In  Parcant  and  N  la  tha  Numbar  of  Polnta  llaad  In  tha  Individual 
Regroaalona 


(deg) 

A 

Tilt  -  0° 
ALfl# 

N 

A 

Tilt  =  20° 

A L„, 

N 

-40 

1.511 

7.2 

31 

1.331 

6.3 

18 

-30 

1.432 

7.2 

36 

1.362 

7.2 

31 

-20 

1.358 

6.9 

36 

1.218 

6.1 

37 

-10 

1.307 

6.8 

37 

1.105 

5.1 

34 

0 

1.262 

5.7 

38 

1.053 

4.5 

32 

10 

1.297 

6.8 

37 

1.074 

4.6 

30 

20 

1.348 

6.9 

36 

1.139 

5.5 

32 

30 

1.429 

7.4 

36 

1.253 

6.0 

27 

40 

1.505 

7.4 

31 

1.291 

5.7 

17 
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which  the  aerosol  would  produce  in  the  single  scatter¬ 
ing  approximation.  This  greatly  simplifies  the  appli¬ 
cation  of  aerosol  models  to  aerosol  analysis  by  satellite 
since  adding  to,  or  in  some  way  changing,  the  aerosol 
model  requires  no  additional  multiple  scattering  com¬ 
putations.  In  fact,  the  only  multiple  scattering  com¬ 
putations  required  for  aerosol  analysis  are  those  for 
determining  I.r,  which  can  be  performed  once  and  for 
all.-9 

The  computations  have  been  explicitly  applied  to 
Band  4  of  the  CZCS,  which,  because  of  its  high  radio- 
metric  sensitivity  and  excellent  calibration,  is  ideal  for 
studying  aerosols  over  the  ocean.  Specifically,  the 
constant  A  in  the  relationship  Lus  =  -  Lr)  is 

given  as  a  function  of  position  along  the  scan  for  four 
typical  orbital-solar  position  scenarios.  For  the  two 
values  of  the  sensor  tilt  angle  examined  (comprising 
~60%  of  the  imagery  acquired  by  CZCS  during  its  first 
year  in  operation),  Las  can  be  retrieved  from  Lt  —  Lr 
with  an  average  error  of  no  more  than  5-7%  except  at 
the  edges  of  the  scan. 

This  work  received  support  from  the  National  Aero¬ 
nautics  and  Space  Administration  under  grant 
NAGVV-273  and  contract  NAS5-28798  and  the  Office 
of  Naval  Research  under  contract  N00014-84-K-0451. 
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